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underwater  distances  by  using  acoustic  modems  and  an  acoustic  ray  tracing  code  to  model 
the  environment  is  proposed  for  real-time  applications. 

To  be  able  to  establish  the  submersible  vehicle’s  position,  a  tracking  algorithm 
relying  on  Kalman  filtering  (KF)  techniques  was  developed  to  fuse  all  available  data.  The 
lack  of  directional  information  from  the  acoustic  modems  employed  produced 
nonlinearities  that  were  treated  using  the  extended  Kalman  filter  (EKF)  and  the  unscented 
Kalman  filter  (UKF). 

The  developed  algorithms  were  initially  tested  using  synthetic  data,  and  the  results 
showed  the  importance  of  a  smoothing  algorithm  to  produce  realistic  trajectories.  This 
analysis  also  suggested  that  faster  convergence  and  better  behavior  in  the  presence  of 
noise  was  achieved  by  the  UKF  approach.  Real  data  collected  during  sea  trials  confirmed 
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I.  INTRODUCTION 


Parts  of  this  chapter  were  previously  published  by  EEEE  [1].* 

There  are  several  definitions  for  navigation.  One  of  our  favorites,  which  was 
translated  into  English,  states  that  “navigation  is  the  science  and  art  of  safely  guiding  a 
vehicle  from  the  starting  point  until  the  final  destination”  [2],  The  simplest  approach  to 
navigating  at  sea  is  based  on  dead  reckoning,  by  which  the  current  vehicle’s  position  is 
estimated  by  integrating  velocity  between  periodic  fixes,  such  as  a  global  positioning 
system  (GPS)  or  similar.  The  inaccuracies  inherent  in  this  process  can  quickly  multiply, 
however,  to  become  significant  errors. 

Advances  in  navigational  aids,  in  particular  satellite  navigation,  have  made  course 
selection  through  the  dead  reckoning  of  human  mariners  obsolete  for  most  purposes; 
although  the  prudent  sea  navigator  still  uses  dead  reckoning  as  a  backup  tool  in  case  of 
failure  of  the  more  advanced  systems. 

In  some  systems,  human  dead  reckoning  has  been  replaced  with  automatic 
(without  human  intervention)  navigational  information  gathered  from  inertial  sensors 
(gyroscopes  and  accelerometers);  such  systems  are  known  as  inertial  navigation  systems 
(INS).  There  are  also  hybrid  systems  that  make  use  of  a  digital  magnetic  compass  and  a 
fusion  of  different  sensors,  or  only  estimates,  to  supply  the  data  necessary  for  navigation. 

In  dead  reckoning  systems,  including  the  INS,  errors  due  to  sensors’  bias  and 
misalignment  accumulate  over  time  and  can  lead  to  unacceptable  position  errors. 
Therefore,  it  is  necessary  to  receive,  regularly,  precise  position  updates  from  external 
reliable  sources,  like  GPS,  to  keep  the  position  errors  at  acceptable  levels. 


*  ©  2016  IEEE.  Personal  use  of  this  material  is  permitted.  Permission  from  IEEE  must  be  obtained  for  all 
other  uses,  in  any  current  or  future  media,  including  reprinting/republishing  this  material  for  advertising  or 
promotional  purposes,  creating  new  collective  works,  for  resale  or  redistribution  to  servers  or  lists,  or  reuse 
of  any  copyrighted  component  of  this  work  in  other  works. 


1 


The  most  advanced  INS  for  underwater  navigation  have  a  margin  of  error  ranging 
from  0.1%  to  1%  of  the  distance  travelled,  without  relying  on  external  positioning 
sources.  This  is  usually  represented  by  what  is  called  the  circular  error  probability,  a 
circle  which  is  bounded  by  50%  probability  of  containing  the  true  position  [3].  From  the 
experiments  conducted  during  this  work,  it  was  found  that  systems  relying  only  on  digital 
magnetic  compasses  and  speed  estimates  may  produce  errors  of  around  20%  of  circular 
error  probability  or  even  higher. 

Successful  surface  navigation  systems  have  been  developed  that  integrate  GPS 
with  inertial  sensors,  but  the  absence  of  GPS  signal  reception  underwater  makes 
navigation  for  manned  and  unmanned  underwater  vehicles  (UUV)  a  more  difficult  task. 
By  surfacing,  the  UUV  can  get  a  position  update  using  its  GPS,  but  this  is  impossible 
(e.g.,  under  ice)  or  undesirable  in  many  circumstances  [4],  In  military  applications, 
surfacing  can  increase  the  chance  of  detection  or  represent  a  delay  in  the  mission  that  is 
unacceptable.  Therefore,  a  technique  to  accurately  update  the  position  of  an  autonomous 
system  while  underwater  is  extremely  desirable. 

Different  approaches  to  establish  an  underwater  positioning  system  have  been 
proposed  over  the  years.  They  include  short  baseline  (SBL),  long  baseline  (LBL),  ultra- 
short  baseline  (USBL),  GPS  intelligent  buoys  (GIB),  and  some  hybrid  systems  based  on 
the  previous  ideas  [5],  [6],  and  [7].  All  of  those  systems  make  use  of  acoustic  travel  time 
measurements  to  estimate  the  distance  and,  in  some  systems,  the  bearing  from  the 
underwater  vehicle  to  reference  points  located  on  the  surface  or  the  sea  floor.  Those 
systems  are  reliable  when  the  depth  being  explored  is  on  the  order  of,  or  larger  than,  the 
horizontal  distances  between  the  assets. 

Estimation  of  horizontal  distances  in  relatively  shallow  water  is  more  complex, 
due  mainly  to  the  multipath  nature  of  the  propagation,  wherein  multiple  arrivals  reach  the 
receiver  at  different  times  by  different  propagation  paths.  Accurately  estimating  the  travel 
time — and,  consequently,  the  distance — of  these  acoustic  signals  under  such  variable 
conditions  is  difficult.  Some  of  the  systems  described  previously  are  able  to  account  for 


2 


the  refraction  of  the  sound  waves  when  the  sound  speed  profile  is  available,  but  none  of 
them  provides  a  solution  to  treat  the  distance  estimation  errors  caused  by  multipath 
signals. 

In  this  dissertation,  the  relatively  shallow  water  environments  wherein  accurately 
estimating  multipath  signals  represents  a  navigational  problem  are  explored.  A  more 
accurate  estimation  for  the  distance,  based  on  acoustic  wave  travel  time  measurements, 
an  acoustic  ray  tracing  code  to  model  the  environment,  and  an  iterative  routine  to  match 
the  measurements  with  synthetic  predictions,  is  proposed.  The  algorithm  is  designed  to 
take  only  a  few  seconds  to  converge  to  a  solution,  making  it  appropriate  for  real-time 
applications  [1]. 

To  obtain  a  reliable  estimate  of  the  vehicle  position,  the  information  for  a  number 
of  sources  and  mathematical  models  must  be  combined  so  effects  due  to  multipath,  sea 
currents,  and  simple  observation  noise  can  be  mitigated.  The  ability  to  establish  the  UUV 
position  via  a  tracking  algorithm  relying  on  Kalman  filtering  (KF)  techniques  is 
developed  to  fuse  all  available  information  (Chapters  II  and  V). 

The  development  of  the  tracking  algorithm  takes  into  account  the  characteristics 
of  the  measurement  system.  This  research  relies  on  battery-powered  digital  signal 
processor  (DSP)-based  acoustic  modems,  which  make  use  of  acoustic  communication 
protocols  to  measure  the  acoustic  wave  travel  time  from  the  UUV  to  the  reference  points. 
The  main  characteristics  of  our  acoustic  modems  are: 

1.  Inability  to  take  simultaneous  measurements  from  different  reference 
points,  and 

2.  No  bearing  (direction)  is  associated  to  the  travel  times  measured  from 
systems  used  in  this  work  (though  such  systems  are  available  for  future 
studies). 

Those  characteristics  produce  a  non-linear  measurement  equation  that  is  treated 
using  two  different  approaches:  a  Taylor  series  linearization  as  in  the  extended  Kalman 
filter  (EKF),  and  a  statistical  linearization  based  on  what  are  called  “sigma  points,”  as  in 
the  unscented  Kalman  filter  (UKF).  Additionally,  practical  aspects  such  as  smoothing  and 
tuning  are  addressed. 
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As  part  of  the  UUV  tracking  model,  the  drift  caused  by  the  sea  current  was 
modeled  as  a  random  walk  and  is  part  of  the  state  of  the  system.  Predictions  for  the  sea 
current  from  different  UUVs  may  be  transmitted,  via  the  reference  point  at  the  surface,  to 
a  command  ship  or  shore-based  command  center  for  further  processing.  In  Chapter  V, 
Section  G,  a  consensus  algorithm  is  presented  to  take  advantage  of  this  information. 

The  developed  algorithms  are  tested  initially  using  synthetic  data  (Chapter  VI), 
and  then  again  using  real  data  collected  during  the  sea  trials  that  took  place  in  Monterey 
Bay  in  August  2015  (Chapter  VII). 

This  dissertation  is  constructed  as  follows.  After  this  Introduction  chapter,  a 
statement  of  the  problem  and  a  model  for  the  UUV  are  presented  in  Chapter  II.  Important 
aspects  of  travel  time  estimation  are  explored  in  Chapter  III,  and  a  technique  to  estimate 
distance  based  on  travel  time  measurements,  employing  a  ray  tracing  code,  is  presented 
in  Chapter  IV.  Two  tracking  algorithms  based  on  the  EKF  and  UKF  are  featured  in 
Chapter  V,  and  simulations  of  the  developed  algorithms  are  discussed  in  Chapter  VI. 
Results  of  the  sea  trials  are  provided  in  Chapter  VII.  A  conclusion  and  recommendations 
for  future  work  are  given  in  Chapter  VIII. 
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II.  PROBLEM  STATEMENT  AND  UUV  MODEL 


As  was  described  in  the  previous  chapter,  underwater  navigation  is  particularly 
challenging  due  to  the  fact  that  a  number  of  navigation  aids,  such  as  GPS,  are  not 
available  underwater.  In  view  of  this,  to  accurately  estimate  a  UUV’s  position  at  any  time 
during  a  mission,  one  may  rely  on  acoustic  communication  between  the  UUV  and  other 
platforms  at  known  locations.  As  addressed  in  later  chapters,  a  combination  of  acoustic 
modems  with  or  without  advanced  features,  such  as  directional  sensitivity,  and  multiple 
surface  vehicles  may  provide  the  required  tracking. 

In  a  typical  mission,  a  UUV  is  part  of  a  network  that  may  consist  of  surface 
platforms  at  known  locations  commonly  provided  by  a  GPS,  bottom  deployed  modems 
(at  known  locations),  and  possibly  other  UUVs.  The  surface  and  bottom  deployed  assets, 
ultimately,  represent  reference  points  from  which  a  UUV,  using  its  acoustic  modem,  can 
estimate  its  distance  and  receive  their  coordinates  (latitude,  longitude,  and  depth),  as 
shown  in  Figure  1. 

Due  to  the  limited  bandwidth  of  the  underwater  channel,  the  typical  acoustic 
modem  is  able  to  communicate  with  only  one  platform  at  a  time.  Thus,  it  is  not  possible 
to  take  measurements  of  the  distance  between  the  UUV  and  two  or  more  reference  points 
at  the  same  time.  Additionally,  we  are  assuming  that  no  bearing  associated  with  the 
distance  measurement  is  available. 
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GPS  satellite 
constellation 


Vehicle  (UUV) 


Figure  1.  Network  of  Reference  Points 


The  main  objective  of  this  dissertation  is  to  improve  the  positioning  accuracy  of 
the  UUV  using  this  network  of  reference  points.  The  goal  is  to  provide  sufficient  tracking 
accuracy  so  that  the  UUV  can  stay  submerged  for  longer  missions  while  still  maintaining 
accurate  track  of  its  position. 

The  approach  to  solve  this  problem  is  to  design  an  algorithm  that  combines  all 
information  available,  which  are  the  travel  time  measurements,  the  acoustic  wave 
propagation  model,  the  UUV  navigation  data  in  terms  of  its  attitude  in  the  water  column 
(heading,  pitch,  depth,  and  speed),  and  the  dynamic  model  of  the  UUV  motion.  All  this 
information  can  be  combined  using  KF  techniques,  in  particular  EKF  and  UKF,  due  to 
the  presence  of  nonlinearities  and  non-Gaussian  disturbances.  The  problem  of  modeling 
is  based  on  the  state  space  representation  for  dynamic  systems  as  can  be  seen  in  the 
next  sections. 

A.  STATE  SPACE  REPRESENTATION  OF  DYNAMIC  SYSTEMS 

A  fundamental  concept  to  describe  the  behavior  of  a  dynamic  system  is  the  state 
of  the  system.  The  concept  of  state  refers  to  a  minimum  set  of  variables  (state  variables) 
that  fully  describe  the  system  and  its  response  to  any  given  set  of  inputs  while  accounting 
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for  errors  due  to  sensor  noise  and  model  uncertainties.  The  set  of  all  possible  values  of 
the  state  variables  is  the  state  space. 

The  state  space  representation  of  a  dynamic  system  is  not  unique  and  usually 
involves  a  time  variable  t  and  three  sets  of  variables: 

State  variables:  xl,x2,  ...,xn 
Input  variables:  uy,u2,  ... ,  um 
-  Output  variables:  z, ,  z2 ,  ...,  z,p. 

For  a  continuous-time  system,  the  description  generally  takes  the  form,  in  vector 
notation,  of  [8] 

x(t)=f(x,u,t)  (2.1) 

z(r)  =  h(x,u,r),  (2.2) 

where  f  (x,u,?)and  h(x,u,?)  are  vector  functions  with  n  and p  components,  respectively. 

Equation  (2.1)  is  called  the  state  equation,  and  Equation  (2.2)  is  called  the  output 
or  measurement  equation.  If  additive  noise  is  present  in  the  state  and  measurement 
equations,  we  have  what  is  called  a  continuous-time  non-linear  stochastic  dynamic 
system,  and  Equations  (2.1)  and  (2.2)  become 

x(t)  =  f  (x,u,t)  +  o(t)  (2.3) 

z(t)  =  h(x,u,t)  +  co(t),  (2.4) 

where  o(t)is  a  n-dimensional  vector  representing  the  input  disturbance  or  process  noise, 

which  is  also  called  plant  noise,  and  ©(?)  is  a  p-dimensional  vector  representing  the 

output  disturbance  or  measurement  noise. 

If  the  system  and  measurements  are  linear,  Equations  (2.3)  and  (2.4)  may  be 
written  in  matrix  notation  as 

(2.5) 

and 
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In  vector  notation  this  becomes 

x(t)  =  A(t)x(t)+B(t)u(t)  (2.7) 

z(f)  =  H(f)x(f)  +  D(f)u(f).  (2.8) 

If  additive  noise  is  present  in  the  state  and  measurement  equations  we  have  what 
is  called  a  continuous-time  linear  stochastic  dynamic  system,  and  Equations  (2.7)  and 
(2.8)  become 

x(t)  =  A(t)x(?)+B(V)u(r)+o(r)  (2.9) 

z  (7)  =  H(t)x(t)  +  D(t)u(t)  +  co(t) .  (2.10) 

For  discrete-time  stochastic  non-linear  dynamic  systems,  the  state  and 
measurement  equations  can  usually  be  represented  by  first  order  difference  equations, 
and  take  the  form  [8] 

x(k  +  l)  =  f  [x(k)]  +  g[u(k)]  +  o(k)  (2.11) 

z(k)  =  h[x(k)]  +  co(k) ,  (2.12) 

where  the  indices  k  and  k+ 1  refer  to  sampling  times  tk=kAt  and  tk+l  =(k  + 1)  At , 

respectively. 

Finally,  for  discrete-time  stochastic  linear  dynamic  systems,  the  state  and 
measurement  equations  usually  take  the  form  [8] 

x(k  +  l)  =  F(k)x(k)+G(k)u(k)  +  v(k)  (2.13) 

z(k)  =  H(k)x(k)  +  D(k)u(k)  +  {o(k).  (2.14) 

B.  UUV  STATE  SPACE  REPRESENTATION 

In  this  section,  a  UUV  state  space  representation  is  presented,  considering  its 
dynamics  and  the  characteristics  of  the  range  measurement  system. 

1.  The  State  Equation 

Assuming  a  UUV  moving  through  the  water  at  depth  d  moving  with  speed  V  , 


we  can  build  a  coordinate  system  where  latitude  and  longitude  are  mapped  to  Cartesian 
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coordinates  ( x  and  y ).  In  this  case,  6  represents  the  heading,  <j>  represents  the  pitch,  and 
d  ( z  axis)  represents  the  depth  of  the  UUV  (see  Figure  2). 


The  y  axis  points  to  east  and  the  x  axis  points  to  north. 

Figure  2.  UUV  Coordinate  System 

The  motion  of  the  UUV  in  the  xy-plane  can  be  represented  by  the  following  set  of 
first  order  differential  equations 

i(?)  =  k(?)cos(0(?))sin(#(?))+c.t  (?)  (2.15) 

y(f)  =  V(f)cos(^(f))cos(^(f))  +  c),  (?) ,  (2.16) 

where  cx  and  cy  represent  the  speed  of  the  sea  current  in  the  x  and  y  directions, 
respectively.  Since  we  assume  the  current  to  be  fairly  constant  in  speed  and  direction  and 
slowly  varying  in  time  and  location,  we  have  chosen  to  model  it  as  a  random  walk  [9]. 
Although  it  is  assumed  to  be  unknown  a  priori,  any  prior  information,  can  be  easily 
incorporated  in  the  proposed  algorithm. 

When  d ,  V ,  0 ,  and  (f)  are  known,  the  dynamic  model  in  its  simplest  form 
becomes 
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In  vector  notation,  Equation  (2.17)  becomes 

x(f)  =  Ax(f)  +  Bu(t) . 


(2.18) 


As  we  are  only  interested  in  knowing  the  state  of  the  system  at  a  discrete  set  of  times 
(sampling  times),  Equation  (2.18)  is  discretized  by  the  following  procedure. 


Performing  a  first  order  Taylor  series  expansion  yields 

x(r+Ar)  =  x(r)  +  x(V)Ar ;  (2.19) 

substituting  (2.18)  into  (2.19)  produces 

x(V  +  Ar)  =  x(t)  +  AAtx(?)  +  BAm(V) ,  (2.20) 

and  rearranging  Equation  (2.20)  results  in 

x(f  +  At)  =  (I  +  AAr)x(r)  +  BAru  (t)  •  (2.21) 

' - - - '  G 

F 

Adding  noise,  we  may  write  the  state  equation  for  the  stochastic  discrete-time 
linear  dynamic  system  as 

x(k  +  l)  =  F(k)x(k)  +  G(k)u(k)  +  o(k).  (2.22) 


In  matrix  notation,  Equation  (2.22)  becomes 
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2.  The  Measurement  Equation 

In  such  a  case  that  the  acoustic  modem  does  not  provide  bearing  information,  we 
have  a  range-only  measurement,  where  each  successful  measurement  represents,  in  the 
xy-plane,  a  circle  of  possible  UUV  positions  (see  Figure  3). 
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Figure  3.  Measurement  Model  in  (a)  yz- Plane  and  in  (b)  xy- Plane 

In  Figure  3,  (xr,yr,zr)  are  the  coordinates  of  the  reference  point,  (x,y,d)  are  the 
coordinates  of  the  UUV,  r  is  the  horizontal  range,  rj  is  the  slant  range,  and  d  is  the  UUV 
depth.  From  Figure  3b,  it  is  noted  that 

(x-xrf  +(y-yrf  =  r2,  (2.24) 

where  is  the  distance  in  the  horizontal  plane. 

Expanding  and  rearranging  Equation  (2.24)  yields 

r2 -x2r -y2r  =x2 -2xxr-2yyr  +  y2 .  (2.25) 

Equation  (2.25)  may  be  written  as 

z(k)  =  h[fc,x(fc)],  (2.26) 

where 

z(k)  =  r(k)2  —xr  (fc)2  —yr  (k)2  is  the  observation,  which  combines  range 
measurements  with  knowledge  of  reference  point  locations; 

h[fc,x(fc)]  =  x(k)2  -2x{k)xr{k)-2y{k^yr  (k)  +  y(k)2  is  a  non-linear 
function. 

With  the  addition  of  noise,  the  measurement  equation  takes  its  final  form  of 

z(k)  =  h[k,x(k)]  +  co(k).  (2.27) 

It  is  important  to  note  that  the  UUV  state  space  representation  is  composed  of  a 

linear  state  equation,  Equation  (2.22),  and  a  non-linear  measurement  equation,  Equation 
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(2.27).  Note,  too,  that  the  measurement  model  needs  only  the  reference  point  coordinates 
and  horizontal  distance  to  the  UUV;  this  means  that  the  reference  point  could  be  a  surface 
asset  able  to  access  a  GPS,  a  bottom  deployed  acoustic  modem  at  a  known  location,  or 
another  UUV  with  a  reliable  position. 


3.  UUV  Model 


The  complete  model  for  the  UUV  dynamics  and  measurements  (state  space 
representation)  is  given  by 


x(k  +  l)  =  F(k)x(k)  +  G(k)u(k)  +  o(k) 
z(&)  =  h[k,x(k)]  +  co(&), 


(2.28) 


where 


The  state  x(k)  is  composed  of  the  UUV  position  and  sea  current,  in  the  xy-plane, 
as  shown  in  Equation  (2.23); 

The  matrices  F  (k)  and  G  (k)  are  deterministic  and  known  for  all  k. 
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The  input  sequence  u(k)is  deterministic  and  known  for  all  k, 
u(£)  =  V  (&)cos(^(A:)) 

z [k)  =  r{k)2  -xr{k^  -yr{k)2  is  the  measurement,  where  r(k)  is  the  horizontal 


sin(#(T)) 
cos  (#(£)) 


distance  from  the  UUV  to  the  reference  point,  the  parameters  xr  (k  )  and  yr(k) 

are  the  coordinates  of  the  reference  point  (longitude  and  latitude  mapped  to 
Cartesian  coordinates); 

h[k,x(k)]  =  x(k)2  -2x(k)xr(k)-2y(k)yr(k)  +  y(k)2  is  a  non-linear  function; 


The  noises  u(k) and  oo(k)  are  assumed  to  be  zero  mean,  white  (i.e., 
uncorrelated)  Gaussian  processes 
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£[°W]  =  E[®W]=0 

«[»(‘)»(i)’]=Q,4i 

E[a(k)a(j)J~\  =  R,StJ 
E  o(k)co(;)T  =0 \/k,j. 


where  Sk  j  = 


1  if  k  =  j 
0  if  k  ^  j 


.  The  UUV  model,  as  developed  in  this  section,  is  the  basis  of 


the  tracking  algorithm  developed  in  Chapter  V. 
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III.  ACOUSTIC  WAVE  TRAVEL  TIME  ESTIMATION 


As  seen  in  Chapter  II,  the  measurement  of  distances  underwater,  from  the  UUV  to 
the  reference  points,  is  an  important  part  of  the  problem  addressed  in  this  dissertation.  In 
this  chapter  and  in  the  next,  the  techniques  used  to  accomplish  this  task  are  discussed. 

Distance  measurement  by  timing  an  echo  is  the  standard  procedure  in  radar  or 
active  sonar  applications.  A  waveform,  with  proper  time  localization  characteristics,  is 
transmitted,  and  the  time  that  it  takes  to  receive  the  echo  of  that  waveform  is  measured 
(two-way  travel  time)  to  estimate  the  distance  from  source  to  target.  The  first  problem  to 
be  solved  is  how  to  estimate  the  travel  time  from  measurements  of  an  acoustic  echo. 

A.  CROSS-CORRELATION  FUNCTION 

It  is  commonly  accepted  that  the  time  of  arrival  of  a  known  pulse  can  be 
determined  by  correlating  the  transmitted  and  received  pulses  [10].  This  is  robust  in  the 
presence  of  noise  and  distortion,  especially  when  the  disturbances  are  Gaussian.  For  two 
jointly  wide-sense  stationary  random  processes,  the  cross-correlation  is  defined  as  an 
ensemble  average: 

R*y(T)=E[x(t)y(t+T)]-  (3>1) 

In  practice,  we  assume  joint  ergodicity  so  that  it  is  computed  by  properly 
averaging  in  time.  Considering  two  ergodic  processes  with  finite  duration  T,  the  time- 
average  cross-correlation  function  may  be  calculated  as 

oo 

Rxy(r)  =  -\x{t)y{t  +  r)dt.  (3.2) 

—oo 

Now  let  us  consider  that  x(t)  =  s(t)  represents  the  transmitted  signal,  defined  in 
the  interval  [0,  T],  and  y(t)  represents  the  received  signal  (single  echo)  defined  in  the 
interval  [rA ,  rA  +  71] ,  where  ta>T. 

For  a  stationary  source  and  target,  the  received  signal  may  be  modeled  as  being  a 
delayed  and  attenuated  version  of  the  transmitted  signal  (plus  noise),  where  y  is  the 


15 


attenuation  factor.  For  simplicity,  dispersion  and  other  propagation  effects  in  the  received 
signal  are  not  being  considered.  Then 

y(r)  =  ys(r-rA  )  +  «(/)•  (3.3) 

The  goal  here  is  to  estimate  the  received  signal’s  time  of  arrival  rA.  From 
Equations  (3.2)  and  (3.3),  the  time-average  cross-correlation  function  may  be  written  as 

Rxy  (r)  =  -  J  5  (0  [yj  (t-rA+T)  +  n(t  +  r)]*  dt .  (3.4) 

—oo 

Considering  that  signal  and  noise  are  uncorrelated  and  evaluating  Equation  (3.4)  at 
t  =  ta,  we  obtain 

■*  00 

K  M  =  Y-  J  K0f  dt  =  Vpavg,s  ■  (3.5) 

—oo 

Equation  (3.5)  shows  that  the  time-average  cross-correlation  function  peaks  at  the 
time  of  arrival  rA  and,  in  this  ideal  scenario,  with  a  value  proportional  to  the  transmitted 

signal’s  average  power  [11].  Therefore,  the  time  when  the  cross-correlation  function 
peaks  is  the  best  estimate  for  the  signal’s  time  of  arrival. 

This  technique  is  usually  applied  in  slowly  varying  environments  where  the 
characteristics  of  the  signal  and  noise  remain  stationary  during  the  finite  observation  time 
T  [12]. 


B.  MATCHED  FILTER 

The  cross-correlation  between  received  and  transmitted  signals  may  be  easily 
implemented  as  a  matched  filter.  The  matched  filter  is  an  optimum  filter  in  the  sense  that 
it  can  maximize  the  signal-to-noise  ratio  (SNR)  at  a  designed  time  instant  t0  after  the 
signal’s  arrival. 

The  impulse  response,  h(t),  of  the  matched  filter  is  defined  by  the  waveform  to 
which  the  filter  is  matched.  Following  [13],  the  impulse  response  of  a  filter  matched  to 
the  signal  s  (?) ,  for  the  case  of  input  white  noise,  is 

hit)  =  Ks*  (ta  —t),  (3.6) 


16 


where  K  is  a  constant  and  t0  is  the  designed  instant  at  which  the  filter  output  will  give  the 
maximum  SNR.  In  applications  that  involve  only  detecting  the  arrival  time  of  a  signal 
immersed  in  noise,  the  choice  for  t0  is  made  only  to  make  the  matched  filter  causal  (if 
necessary);  that  is,  /z(t)  must  be  zero  for  t  <0.  It  implies  that  the  parameter  t„  shall  be 

greater  than  or  equal  to  the  duration  of  the  transmitted  signal  T  [14].  If  causality  is  not  a 
concern  in  the  matched  filter  implementation,  ta  can  be  set  to  zero,  as  in  the  example 
shown  in  Figure  4. 

For  completeness,  the  frequency  response  of  the  matched  filter  is 

H(f)  =  KS*(f)e~J2*ft°.  (3.7) 

The  matched  filter  output  y(t)  may  be  represented  by  the  convolution  integral  of 
the  received  signal  (as  in  Equation  (3.3)  and  the  filter  impulse  response  [14] 

oo 

y(t)  =  |  [ys(u - ta)  +  n(u.y^s* (u -  t  +  ta)du ,  (3.8) 

—00 

where  K  in  Equation  (3.6)  was  set  to  1.  Given  that  the  signal  and  noise  are  uncorrelated, 
Equation  (3.8)  reduces  to 

oo 

y(t)=  |  y s(u-TA)s*(u-t  +  t0)du .  (3.9) 

—00 

Equation  (3.9)  may  be  recognized  as  the  cross-correlation  of  the  received  signal  with  the 
transmitted  signal  lagged  by  ta-t.  The  peak  in  the  matched  filter  output  will  occur  when 
u-rA=u-t+t0,  which  means  that  tpeak  =  rA  +  tQ  because 

00  00 

y(TA+0=  j  ys(u-TA)s*(u-TA)du  =  yj  \s(t)fdt.  (3.10) 

—00  —00 

Comparing  Equation  (3.5)  with  Equation  (3.10)  we  can  see  that  both  give  similar  results. 
At  this  point  one  may  conclude  that  the  matched  filter  and  the  cross-correlation  function 
are  equivalent  ways  to  estimate  a  signal’s  time  of  arrival. 
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Transmitted  Signal 


Received  Signal 


Received  Signal+noise 

10  r 
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time  [s] 


Matched  Filter  Output 


Figure  4.  Matched  Filter  Used  to  Detect  Signal’s  Time  of  Arrival, 

Designed  with  ta= 0 

Figure  4  represents  the  output  of  a  matched  filter  designed  with  t0= 0.  In  this 
example,  a  hyperbolic  frequency  modulated  pulse  is  transmitted  (sweeping  from  200  to 
2200  Hz,  7=  50  ms,  7F>=100)  and  a  signal  immersed  in  noise,  SNR=3dB,  is  received  (left 
panel,  first  and  third  plots).  In  an  ideal  no-noise  scenario,  the  received  signal  arrives  at 
t=0.ls  (left  panel,  second  plot).  The  matched  filter  peaks  at  t=0.ls  for  both  signals  (no 
noise  and  signal  plus  noise),  indicating  a  reception  at  this  time  (right  panel). 

1.  The  Doppler  Effect 

When  there  is  relative  motion  between  source  and  target,  the  performance  of  the 
matched  filter  degrades.  Referring  to  Figure  5,  consider  the  simple  case  when  a  source  in 
motion  at  constant  velocity  is  transmitting  a  continuous  wave  (CW)  pulse  during  7 
seconds,  and  receiving  an  echo  from  a  target  in  motion.  The  relative  motion  causes  a  shift 
in  the  frequency  of  the  received  signal.  This  effect  is  called  the  Doppler  effect. 
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Vr=  I  V,||  -  VrMl 


Source 


Figure  5.  Relative  Radial  Speed  between  Source  and  Target 


Defining  the  Doppler  shift  as  fd=  fc—  fR  where  fc  is  the  frequency  of  the 
transmitted  signal,  fR  is  the  frequency  of  the  received  signal,  and  assuming  that  the 
propagation  speed  in  the  medium  (c)  is  much  greather  than  the  relative  radial  speed 
between  source  and  target  (vr)  yields  [15] 


(3.11) 


c 


This  frequency  shift  may  affect  the  system’s  ability  to  detect  the  incoming  signal,  which 
is  why  it  is  necessary  to  take  it  into  account  when  designing  the  matched  filter.  Especially 
in  underwater  applications,  while  the  relative  motion  might  be  small  compared  to  the 
medium  propagation  speed,  the  carrier  frequency  (for  modulated  waveforms)  might  be 
sufficiently  high  so  that  the  frequency  shift  becomes  relevant. 

2.  Received  Signal  Model 

To  understand  the  effect  of  Doppler  shift  in  the  received  signal,  let  us  consider  a 
target  moving  with  constant  radial  speed  vr  away  from  a  static  source.  At  time  t  the 

target  will  be  located  at  range  r(f) ,  given  by 


r(t)  =  ro+vrt. 


(3.12) 
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where  at  t=0  the  target  will  be  located  at  range  r  .  Assuming  that  the  source  will  start 

transmitting  x(t)  at  t  =  0,  the  received  signal  (echo)  at  time  t  is 

x'(t)  =  x[t-r(t)] .  (3.13) 

As  the  target  is  moving  while  the  signal  is  being  reflected,  the  two-way  travel  time  r(t) 

in  Equation  (3.13)  is  considered  a  function  of  time  [16].  For  simplicity,  no  attenuation  or 
other  propagation  effects  are  being  considered. 


Assuming  that  the  propagation  speed  in  the  medium  c  is  constant,  we  can 
calculate  r(t)  as  [17] 

r  (O' 


r(t)  =  -r 


t-- 


(3.14) 


Substituting  Equation  (3.12)  into  Equation  (3.14)  yields 


/  ,  2 
r(,)  =  - 


r + vt  —  - 


VT 


(O' 


2r„+2lv 

c  +  v,. 


(3.15) 


Then  by  substituting  Equation  (3.15)  into  Equation  (3.13),  we  obtain 


:'«  =  ■ 


t- 


2  ro  +  2 vrt 
c  +  v. 


=  x 


J 


1- 


2v. 


y 


c  +  v.. 


t- 


VV  "  '  vrJ 

From  Equation  (3.16),  the  received  signal  may  be  written  as 


2  r 


c  +  v. 


r  J 


:'(t)  =  x(at-T0), 


(3.16) 


(3.17) 


2v  2  r 

where  for  c » vr ,  a«l - -  and  to  «  — - . 

c  c 

As  represented  in  Equation  (3.17),  the  Doppler  effect  manifests  by  compressing 
(or  expanding)  the  received  signal,  and  this  is  mathematically  described  by  a  scale  factor 
( a )  in  time.  Even  if  the  transmitted  signal  can  be  considered  narrowband,  the  effect  in 
the  received  signal  is  a  translation  in  the  frequency  of  the  carrier  (Doppler  shift)  and  an 
increase  or  decerease  in  the  signal  bandwidth. 

Due  to  the  frequency  mismatch  between  transmitted  and  received  signals,  a  filter 
matched  to  v(t)  will  no  longer  be  matched  to  x'(f)  in  Equation  (3.17).  This  mismatch 


may  cause  a  reduction  in  the  peak  at  the  matched  filter  output,  and  depending  on  the 
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severity  of  the  mismatch,  this  peak  can  be  reduced  to  levels  that  do  not  permit  target 
detection. 

If  the  relative  radial  speed,  vr ,  is  known  or  the  Doppler  shift  can  be  measured,  a 
filter  matched  to  x'(t )  may  be  constructed  as 

h{t)  =  x[a(t0-t)\.  (3.18) 

In  cases  where  the  relative  radial  speed  is  not  known  in  advance,  the  classical 
approach  is  to  have  a  multi-hypothesis  detector,  also  called  a  “bank”  of  matched  filters, 
tuned  at  different  radial  speeds  (or  Doppler  shifts),  where  the  filter  that  gives  the  highest 
output  is  selected  (see  Figure  6).  This  process  is  optimal  under  restricted  circumstances, 
computationally  intense,  and  time  consuming. 


Select  Max. 
Amplitude 


Figure  6.  Multi-hypothesis  Detector  or  “Bank”  of  Matched  Filters 

The  use  of  such  detectors  can  introduce  serious  limitations  in  some  real-time 
underwater  acoustics  applications  that  take  advantage  of  physically  small,  battery 
powered  DSP-based  modems  [18].  In  these  applications,  the  use  of  a  signal  with  a  certain 
tolerance  for  the  Doppler  effect  is  desirable  to  avoid  the  use  of  complicated  and  time 
consuming  detectors. 
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C.  COMMON  WAVEFORMS  IN  UNDERWATER  ACOUSTICS 

The  selection  of  the  waveform  is  a  complex  and  important  task,  playing  an 
important  role  in  the  overall  performance  of  distance  measurements  of  systems 
underwater.  The  most  common  waveforms  used  for  this  task  are  the  linear  frequency 
modulated  (LFM)  pulse  and  the  hyperbolic  frequency  modulated  (HFM)  pulse. 

The  popularity  of  these  signals  comes  from  their  pulse  compression 
characteristics.  Pulse  compression  allows  achieving  transmitted  power  of  a  relatively 
long  pulse  while  obtaining  the  range  resolution  corresponding  to  a  short  CW  pulse  [19]. 

These  two  classes  of  signals,  are  defined  in  subsections  (a)  and  (b). 
a)  LFM  pulse  [20] 

x(t}  =  a(t^cos[27rfct  +  Dpt2^  for  |t|<T/2,  (3.19) 

where  a(t)is  a  real  amplitude  modulating  function,  Dp=±^ B  is  the  phase  deviation 

constant  in  rad/sec2,  B  is  the  swept  bandwidth  in  Hz,  and  T  is  the  pulse  duration  in 
seconds.  The  instantaneous  frequency  is  defined  as 

In  complex  notation,  Equation  (3.19)  then  becomes 

x(t)  =  a{t)ei(2*U+D/\  (3.21) 


) 

(3.20) 
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This  LFM  pulse  is  defined  in  the  interval  0  <  t  <  50ms  with  B  =  200Hz  (77?=  10). 

Figure  7.  Rectangular-Envelope  (a)  LFM  Pulse  and  (b)  Instantaneous  Frequency 


b)  HFM  pulse  [21] 

x(t )  =a(t)  cos 


—  ln(l  +  Kfat) 

K 


for  \t\  <  r/2 , 


(3.22) 


where  a(t)is  a  real  amplitude  modulating  function,  k  = 


f o  fend 


>  fo  aild  fend  ai'e’ 


ffofend 

respectively,  the  starting  and  the  ending  frequency  in  Hz,  and  T  is  the  pulse  duration  in 
seconds.  The  HFM  swept  bandwidth  is  defined  as  B  =  +(fo  -  fend  )  (+  for  up  sweep  and  - 
for  down  sweep).  The  instantaneous  frequency  is 


j,  /  \  1  d 

dt 

 fo 


— ln(l  +  Kf„t) 

K 


(3.23) 


1  +  Kfot 

In  complex  notation,  Equation  (3.22)  then  becomes 


v(t)  =  a(t)e 


j—\n(l+Kf0t) 


(3.24) 
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This  HFM  pulse  is  defined  in  the  interval  0  <  t  <  50ms  with  B  =  200Hz  (T6=10). 

Figure  8.  Rectangular-Envelope  (a)  HFM  Pulse  and  (b)  Instantaneous  Frequency 


The  behavior  of  LFM  pulses  in  the  presence  of  Doppler  will  depend  on  the  signal 
time-bandwidth  (TB)  product,  source-target  relative  radial  speed  (vr),  and  the 

characteristic  medium  propagation  speed  ( c ).  In  a  stationary  source  and  target  scenario, 
the  nominal  half- width  of  the  main  lobe  in  the  matched  filter  output  is  1/B  [17].  The 
Doppler  will  cause  a  widening  in  the  main  lobe  that  is  accompanied  by  a  reduction  in 
amplitude  (due  to  conservation  of  energy)  of  the  matched  filter  output’s  peak.  Depending 
on  how  severe  the  widening  effect  is,  it  will  be  necessary  to  use  a  “bank”  of  matched 
filters  to  permit  target  detection. 

The  matched  filter  output  is  shown  in  Figure  9  for  three  scenarios  of  source-target 
relative  motion:  stationary,  relative  radial  speed  of  7.5  knots,  and  relative  radial  speed  of 
15  knots.  The  transmitted  waveform  is  a  rectangular-envelope  LFM  pulse  of  50  ms 
duration,  sweeping  from  9  kHz  to  14  kHz  (77?=250).  As  shown  in  Figure  9,  when  there  is 
no  motion,  the  pulse  arrives  at  0.1s. 
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The  characteristics  of  the  LFM  pulse  are:  T=50ms,  B=5  kHz  sweeping  from  9  kHz  to  14 
kHz  (TB=250). 

Figure  9.  Matched  Filter  Output  Having  as  Input  a  Rectangular-Envelope  LFM  Pulse 

for  Different  Source-Target  Relative  Radial  Speeds 


Note  that  the  width  of  the  main  lobe  increases  with  the  increase  of  the  radial 
speed,  and  the  peak  amplitude  decreases.  Note,  too,  that  there  is  a  temporal  offset  in  the 
matched  filter’s  output  peak  when  relative  motion  between  source  and  target  is  present. 
Figure  9  illustrates  a  negative  Doppler  shift  (source  and  target  are  receding)  so  the 
relative  speed  increases  the  round  trip  time  delay  estimated  at  the  receiver. 

As  demonstrated  in  the  Appendix,  the  widening  in  the  main  lobe  of  the  matched 
filter  output  is  given  by  the  factor  l  +  (4 vr/c)TB .  If  (4 vr/c)TB  can  be  considered  much 

smaller  than  unity,  the  widening  and,  consequently,  the  amplitude  reduction  in  the 
matched  filter  output  peak  may  be  neglected,  which  implies 

TB  «  —  .  (3.25) 

4  vr 
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Equation  (3.25)  states  that  if  the  signal’s  time  bandwidth  product  is  small  in 
comparison  with  the  ratio  of  medium  sound  speed  over  the  maximum  source-target 
relative  radial  speed,  the  widening  effect  in  the  matched  filter  output  may  be  neglected.  In 
radar  and  laser  applications,  due  to  the  large  difference  between  c  andvr,  LFM  pulses 

can  be  used  without  issues,  even  for  relatively  large  TB.  However,  in  underwater 
acoustics,  due  to  the  relative  lower  sound  speed  propagation,  the  use  of  this  signal  can 
introduce  limitations  in  the  system’s  ability  to  detect  a  target. 

For  HFM  pulses  (as  demonstrated  in  the  Appendix),  the  main  lobe  in  the  matched 
filter  output,  in  the  presence  of  Doppler,  is  slightly  different  from  1/B  and  independent  of 
(4 vJc)TB .  Thus,  the  widening  and,  consequently,  the  reduction  in  the  matched  filter 

output’s  peak  will  not  be  as  severe  as  in  the  LFM  case,  and  may  permit  target  detection 
without  the  use  of  a  multi-hypothesis  detector. 

The  matched  filter  output  is  shown  in  Figure  10  for  three  scenarios  of  source- 
target  relative  motion:  stationary,  relative  radial  speed  of  7.5  knots,  and  relative  radial 
speed  of  15  knots.  The  transmitted  waveform  is  a  rectangular-envelope  HFM  pulse  of 
50  ms  duration,  sweeping  from  9  kHz  to  14  kHz  (7B=250),  and  a  single  echo  is 
considered  to  be  received  at  t=0.1s.  Note  that  the  width  of  the  main  lobe  barely  changes 
with  the  increase  of  the  radial  speed,  and  the  peak  amplitude  slightly  decreases.  Note,  too, 
that  as  in  the  LFM  case,  there  is  a  temporal  offset  in  the  matched  filter’s  output  peak 
when  relative  motion  between  source  and  target  is  present. 
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The  characteristics  of  the  HFM  pulse  are:  T=50ms,  B=5  kHz  sweeping  from  9  kHz  to  14 
kHz  (TB=250). 

Figure  10.  Matched  Filter  Output  Having  as  Input  a  Rectangular-Envelope 
HFM  Pulse  for  Different  Source-Target  Relative  Radial  Speeds 


As  seen  in  Figures  9  and  10,  with  the  use  of  both  waveforms,  the  relative  speed 
between  source  and  target  causes  a  temporal  offset  in  the  matched  filter  output’s  peak. 
Compensation  for  this  temporal  offset  should  be  applied  to  obtain  a  precise  travel  time 
measurement. 

The  use  of  HFM  pulses  permits  a  direct  compensation  for  this  time  offset  by  using 
a  closed  form  expression  (Appendix,  Equation  A.61)  when  the  Doppler  shift  is  known. 
Unfortunately,  there  is  no  closed  form  for  LFM  pulses. 

Due  to  these  effects,  most  ranging  applications  in  underwater  acoustics  use  HFM 
pulses  and,  as  will  be  seen  in  the  next  chapter,  HFM  pulses  are  employed  in  the  acoustic 
modems  used  in  this  work. 
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IV.  RANGE  ESTIMATION 


Parts  of  this  chapter  were  previously  published  by  IEEE  [1].* 

As  seen  in  Chapter  III,  the  standard  technique  to  estimate  the  distance  to  a  target 
in  radar  or  active  sonar  applications  begins  with  the  determination  of  the  travel  time.  It 
can  be  done  by  cross-correlating  the  received  signal  (the  backscattered  echo)  with  a 
replica  of  the  transmitted  signal  (or  match  filtering  the  received  signal).  The  two-way 
travel  time  can  then  be  estimated  from  the  time  lag  where  the  peak  in  the  cross¬ 
correlation  occurs. 

In  this  chapter,  a  technique  to  estimate  the  distance  using  the  measured  travel  time 
and  a  ray  tracing  code  to  model  the  acoustic  channel  is  addressed.  To  start  this 
discussion,  the  next  section  provides  a  basic  description  of  the  ranging  routine  used  by 
the  Teledyne-Benthos  acoustic  modems. 

A.  ACOUSTIC  MODEM  RANGING  ROUTINE 

In  this  dissertation,  Teledyne-Benthos  ATM-900  series  acoustic  modems  have 
been  used.  The  modems  installed  in  the  different  assets  have  a  built-in  ranging  routine 
that  makes  use  of  HFM  pulses  to  estimate  the  two-way  travel  time  between  modems.  The 
HFM  pulse  used  by  the  acoustic  modems  has  a  50  ms  duration  and  5  kHz  bandwidth 
(sweeping  from  9  to  14  kHz).  The  basics  of  the  Teledyne-Benthos  built-in  ranging 
routine  can  be  described  in  the  following  manner. 

Consider  a  situation  where  two  modems  are  separated  by  some  distance  R  and 
modem- 1  requests  a  range  to  modem-2,  as  indicated  in  Figure  11. 


_  ©  2016  IEEE.  Personal  use  of  this  material  is  permitted.  Permission  from  IEEE  must  be  obtained  for  all 
other  uses,  in  any  current  or  future  media,  including  reprinting/republishing  this  material  for  advertising  or 
promotional  purposes,  creating  new  collective  works,  for  resale  or  redistribution  to  servers  or  lists,  or  reuse 
of  any  copyrighted  component  of  this  work  in  other  works. 
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Modem-2 


Modem- 1 


R 


Modem-2 


Modem- 1 


Figure  11.  Teledyne-Benthos  Ranging  Routine.  Adapted  from  [22], 

Initially,  modem- 1  transmits  via  MFSK  (multiple  frequency  shift  keying)  a  utility 
package  to  modem-2  preceded  by  an  HFM  pulse,  recording  the  time  when  the  routine 
starts  (tB).  The  HFM  pulse  can  be  considered  as  a  probe  signal  that,  after  matched 

filtering,  provides  an  estimate  of  the  channel  impulse  response  (Section  B.2.b).  With  the 
help  of  a  matched  filter  in  modem-2,  the  incoming  signal  is  detected  and  an  estimate  for 
the  time  of  arrival  is  made  based  on  the  highest  peak  in  the  matched  filter  output. 
Modem-2  then  replies  to  the  range  request  after  a  known  time  delay  td  ,  sending  an 

MFSK  message  (“response”)  containing  the  time  delay  (and  other  information  unrelated 
to  range  estimation)  preceded  by  an  HFM  pulse.  Modem- 1  now  estimates  the  time  of 
arrival  (tE )  from  the  highest  peak  in  the  matched  filter  output.  From  the  time  difference 

between  tE  and  tB  (minus  zD  ),  the  two-way  travel  time  ( r2way )  is  estimated  by  modem- 1 
according  [22],  [23] 

*2 way  ~  i?E  —  )  _  •  (4- 1 ) 

At  this  point,  modem- 1  provides  a  rough  estimate  of  the  slant  range  R  by 
multiplying  the  one-way  travel  time  tm  =  r2way  j 2 ,  by  the  characteristic  sound  speed  of  the 

medium.  A  more  accurate  technique  to  estimate  the  distance  using  a  ray  tracing  code  to 
model  the  environment  is  addressed  in  Section  B.2. 

For  completeness,  it  is  worth  noting  the  capability  of  the  modems  to  compensate 
for  the  Doppler  effect.  The  ability  of  HFM  pulses  to  enhance  target  detection  in  the 
presence  of  Doppler  shift  was  discussed  in  Chapter  III  and  is  demonstrated  in  the 

30 


Appendix.  In  addition,  a  closed  expression  for  the  temporal  offset  in  the  matched  filter 
output  due  to  source-target  relative  motion  is  provided. 

By  measuring  the  Doppler  shift,  the  temporal  offset  in  the  matched  filter  output 
can  be  compensated  using  a  variation  of  Equation  (A.61)  of  the  Appendix.  Teledyne- 
Benthos  acoustic  modems  have  a  patented  algorithm  that  makes  use  of  single  frequency 
tonals  to  measure  the  Doppler  shift  and  then  compensate  for  the  temporal  offset  in  the 
matched  filter  output  [18]  without  relying  on  a  new  matched  filter. 

B.  DISTANCE  MEASUREMENT 

In  an  environment  where  the  speed  of  propagation  is  approximately  constant,  the 
distance  between  source  and  target  may  be  estimated  by  a  simple  multiplication  between 
the  characteristic  propagation  speed  of  the  medium  and  the  estimated  one-way  travel 
time.  This  is  generally  not  the  case  in  underwater  acoustics,  however,  where  the  sound 
speed  varies  spatially.  Even  in  a  benign,  range-independent  environment,  the  sound  speed 
is  generally  observed  to  vary  with  depth  (see  Figure  12). 


Sound  Speed  (m/s) 

Figure  12.  Sound  Speed  Profile  Measured  in  Monterey  Bay  on  August  12,  2015 
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Another  important  characteristic  in  underwater  acoustics  is  the  multipath  nature  of 
the  propagation.  In  a  multipath  environment,  multiple  arrivals  reach  the  receiving 
transducer  at  different  times  by  different  propagation  paths,  causing  several  peaks  in  the 
matched  filter  output.  A  common  approach  is  to  estimate  the  travel  time  by  the  time  of 
the  highest  peak,  which  considers  that  the  highest  peak  represents  the  arrival  with  higher 
energy  and,  therefore,  the  direct  path  between  source  and  target  (see  Section  A). 

Occasionally,  in  a  multipath  environment  where  the  horizontal  distance  between 
source  and  target  is  much  larger  than  the  depth  of  the  water,  the  highest  peak  in  the 
matched  filter  output  (representing  the  arrival  with  higher  energy)  does  not  represent  a 
wave  traveling  along  the  direct  path.  For  example,  the  highest  peak  could  be  due  to  an 
arrival  taking  the  bottom/surface  bounce  path,  as  depicted  in  Figure  13.  An  attempt  to  use 
the  measured  travel  time  to  estimate  a  distance  between  source  and  target,  assuming  that 
it  is  representative  of  the  direct  path,  will  cause  errors  in  the  final  distance  estimation. 


Range  (m)  time  (s) 

In  ray  tracing,  eigenrays  are  defined  as  rays  that  connect  the  source  to  a  particular 
receiver  location.  Note  that,  in  this  example,  the  first  group  of  arrivals  has  smaller 
amplitude  than  the  second  group. 


Figure  13.  Example  Calculation  of  (a)  Eigenray  Paths  and 

(b)  Eigenray  Times  of  Arrival  for  Omnidirectional 
Source  and  Receiver  Using  BELLHOP  [24] 
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In  this  dissertation,  an  approach  is  developed  that  utilizes  an  efficient  ray  tracing 
code  [24],  in  situ  sound  speed  measurements,  sea  floor  modeling,  and  the  measured 
acoustic  wave  travel  time  to  estimate  the  distance  between  the  UUV  and  the  reference 
points.  The  ray  tracing  code  can  then  distinguish  cases  in  which  the  strongest  arrival  is 
not  representative  of  the  direct  path  (or  first  arrival).  Additionally,  with  information  of 
the  sound  speed  profile,  the  refraction  of  the  sound  waves  is  accounted  for  in  the  distance 
estimation. 


1.  Acoustic  Ray  Paths 

Initial  attempts  at  modeling  sound  propagation  underwater  were  motivated  by 
problems  in  predicting  sonar  performance  during  World  War  II  in  support  of  anti¬ 
submarine  warfare  operations.  These  first  models  used  ray  tracing  techniques  derived 
from  the  wave  equation  [25], 

From  a  historical  point  of  view,  the  behavior  of  ray  paths  was  studied  long  before 
acoustic  ray  theory  was  mathematically  formalized.  Ray  theory  originally  emerged  from 
optics,  where  it  was  used  to  understand  propagation  of  light  even  before  the  Maxwell 
equations  were  known.  In  fact,  ray  paths  were  originally  studied  by  Euclid  (around  300 
BC)  [26]  and  complemented  by  the  Arabic  scientist  Alhazen  in  the  tenth  century,  while 
Snell’s  law  only  dates  back  to  1626  [27]. 

Nowadays,  the  research  community  has  lost,  in  part,  some  interest  in  ray  tracing 
codes  mainly  due  to  their  inherent  high  frequency  approximation,  which  in  certain  cases 
can  lead  to  coarse  accuracy  in  the  results.  In  high  frequency  applications,  however,  ray 
tracing  algorithms  are  very  popular  due  to  their  accuracy  and  high  processing  speeds.  A 
brief  discussion  about  the  high  frequency  approximation  can  be  found  in  Section  B.l.b. 

a.  The  Acoustic  Wave  Equation 

The  linearized  wave  equation  can  be  derived  from  the  equation  of  continuity,  the 
Euler  equation,  and  the  equation  of  state  [28].  In  the  absence  of  density  discontinuities,  it 
takes  the  form 
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(4.2) 


vV 


1  d2p 

-flit2 


=  0. 


Density  discontinuities  at  interfaces  can  be  treated  separately  within  the  context  of 
reflection  and  transmission  coefficients  at  boundaries. 


To  obtain  a  time-independent  wave  equation  (Helmholtz  equation),  a  time 
harmonic  solution  for  the  acoustic  pressure  p  is  assumed  at  the  form 

p  =  p(r)eiM,  (4.3) 

where  the  position  vector  r  extends  from  the  reference  point  to  a  certain  point  in  the 
field.  Substituting  Equation  (4.3)  into  Equation  (4.2),  the  wave  equation  reduces  to  the 
Helmholtz  equation 

V2/?(r)  +  &(r)2  p(r)  =  0,  (4.4) 

where  k(r)  =  co/c(r)  An  this  form,  all  of  the  spatially  varying  environmental  factors  that 
affect  the  propagation  are  incorporated  into  k  (r  ) . 


b.  Mathematical  Derivation 

Assuming  a  solution  in  the  format  of  the  product  of  an  amplitude  function  A(r) 
and  a  phase  function  <z»r(r) ,  we  may  define 


p(r)  =  A(r)e7<w4r)  (4  5) 

Substituting  this  into  the  wave  equation  and  equating  the  real  and  imaginary  parts,  yields 


1  '  c(rf  A(f )<r  ' 


and 


2  VA(r)  - Vr(r)  + A(r)  V2r(r)  =  0 . 
Taking  the  asymptotic  limit  co  — >  co  in  Equation  (4.6),  we  obtain 

Mf)|2=-4t 

cy) 


(4.6) 


(4.7) 


(4.8) 
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Equation  (4.8)  is  referred  as  the  Eikonal  equation  and  describes  the  evolution  of 
the  phase  front  in  space  or  the  ray  paths  (see  Figure  14).  Equation  (4.7)  is  called  the 
Transport  equation  and  describes  the  evolution  of  the  amplitude  along  the  ray  paths. 


Figure  14.  Rays  and  Phase  Fronts.  Adapted  from  [26]. 


Instead  of  assuming  the  asymptotic  limit  ru— » oo  in  Equation  (4.6),  a  different 

V2A(F )  1 

argument  can  be  made.  For  frequencies  such  that  — .  .  ;  <sc - T ,  the  Eikonal 

A(r)co  C(F)2 

equation  is  achieved. 


Other  authors  propose  a  similar  argument  as  in  [28],  requiring  that 


v2a(f) 

A(F) 


«;  —  2  .  According  to  this  argument,  the  sufficient  condition  to  be  satisfied  is 

c(f) 


that  the  amplitude  of  the  wave  and  the  speed  of  sound  must  not  change  significantly  over 
distances  comparable  to  a  wavelength.  The  result  of  all  of  the  previous  arguments  is 
basically  the  same,  and  at  this  point  we  can  just  say  that  a  high  frequency  approximation 
is  used  to  obtain  the  Eikonal  equation. 


c.  Solving  the  Eikonal  Equation 

The  Eikonal  equation  is  a  first-order,  nonlinear,  partial  differential  equation  that  is 
usually  solved  by  the  method  of  characteristics  [29].  The  method  of  characteristics 
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essentially  involves  recognizing  that  V  r  is  everywhere  perpendicular  to  the  phase  front 
[30].  This  defines  the  local  direction  of  the  ray  trajectories,  which  are  everywhere 
perpendicular  to  the  acoustics  phase  fronts. 

If  we  introduce  a  family  of  curves  that  describes  the  evolution  of  the  phase  front 
by  defining  them  as  perpendicular  to  the  phase  front  everywhere,  then  the  ray  paths  are 
points  in  space,  r(s),  defined  by 

-  =  cVr,  (4.9) 

ds 

where  s  is  the  distance  along  the  trajectory  or  the  arc  length  along  the  ray.  The  factor  c  is 
included  in  Equation  (4.9)  such  that  the  tangent  vector  along  the  ray  has  unit  length. 


It  can  be  shown  [26]  that  solutions  to  the  ray  paths  in  two  dimensions  satisfy  a  set 
of  coupled,  first-order  differential  equations  of  the  form 


dr  dz  n 

—  =  ac ,  —  =  Be , 
ds  ds 


(4.10) 


da  _  1  dc  d/3  _  1  dc 
ds  c2  dr  ’  ds  c2  dz 


(4.11) 


that  can  be  numerically  solved  by  a  variety  of  methods.  Considering  a  trajectory  in  two 
dimensions  (r,z),  we  require  as  an  initial  condition  that  the  rays  start  at  the  source 

position  ( r  ,  z„ )  and  have  a  defined  launch  angle,  0o  (see  Figure  15),  where 


tan{0o)  =  ^=di^  =  ^p-. 
dr  dV/ds  C0SW 

Substituting  Equation  (4.10)  into  Equation  (4.12),  we  obtain 

sin^=js 


(4.12) 


(4.13) 


COHA,)  =  —  =  aroC(KrZo) 


1  dr 
C{ro’Zo)  ds 


(4.14) 


which  provides  the  initial  conditions  for  the  parameters  a  and  (3  in  terms  of  the  launch 
angle  and  the  speed  of  the  sound  at  the  source  position. 
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Figure  15.  Rays  2D  Initial  Conditions 


The  ray  paths  characterized  by  Equations  (4.10)  and  (4.11),  along  with  the  initial 
conditions,  describe  the  direction  in  which  the  phase  front  evolves  spatially.  But  they  do 
not  describe  the  amplitude  of  the  pressure  field.  To  obtain  the  complete  description  of  the 
pressure  field,  it  is  necessary  to  associate  a  phase  and  amplitude  for  each  ray. 

The  phase  along  a  ray  path  can  be  determined  from  the  Eikonal  equation, 


Vr- Vr  =  4r  • 


(4.15) 


Substituting  Equation  (4.9)  into  Equation  (4.15),  we  obtain 


c  ds  c1 


and  simplifying  yields 


dr  _  1 
ds  c 


'(Hw)ds' ' 


(4.16) 


(4.17) 


where  r(.s)  is  the  travel  time  along  the  ray  and  the  phase  is  cot(s).  In  a  multipath 
environment,  7(5)  is  important  to  estimate  the  time  which  different  arrivals  reach  the 
receiver,  as  depicted  in  Figure  13b. 
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d.  Solving  the  Transport  Equation 


The  amplitude  along  a  ray  path  is  obtained  by  solving  the  Transport  equation, 
which  may  be  written  as 

V-(a(F)2Vt(F))  =  0.  (4.18) 

The  solution  to  Equation  (4.18)  can  be  developed  in  terms  of  ray-tube  areas  and 
conservation  of  energy  flux.  Consider  a  ray  propagating  in  range  from  ro  to  r,  as 
depicted  in  Figure  16,  a  ray  tube  can  be  constructed  from  nearby  rays  passing  through  the 
small  area  S0 .  When  the  ray  tube  reaches  r ,  its  cross  sectional  area  will  be  S . 

Integrating  Equation  (4.18)  over  the  volume  of  the  ray  tube  segment,  applying  Gauss’ 
theorem  and  recognizing  that  the  rays  are  normal  to  the  phase  fronts  and  therefore  vanish 
on  the  sides  of  the  ray  tube  yields  [26],  [31] 


and 


A(r0f 

pc  P0c0 


(4.19) 


A(r)  =  A(rc) 


(4.20) 


Equation  (4.20)  implies  a  rising  or  falling  amplitude  as  the  ray  tube  shrinks  and  expands. 


Figure  16.  Geometry  of  a  Ray  Tube 
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As  shown  in  [30],  the  amplitude  of  the  field  within  a  ray  “tube”  can  be  expressed 


as 


A(r)  =  A(r) 

V  r  poca  8(j)  sin  9  8z 


U(oJ 


r„  pc  8(j)o  sin  0o  ro80o 


r  poco  <50  sin#  8z 


In  terms  of  transmission  loss  ( TL ), 

TL  =  -  lOlog 


f 

|2  'N 

k!_ 

P2 

V 

0  J 

or 


TL(r,z)  =  lOlog 


f  ^ 

r 


\roJ 


-lOlog 


f  Pc  ^ 
P  c 

r o  o 


lOlog 


dB  re  r  , 


^  <50,  sin  6o  x 
8(f)  sin  9 


+  lOlog 


8z(r,z) 

V  roS6o  J 


(4.21) 


(4.22) 


(4.23) 


The  first  term  in  Equation  (4.23)  represents  cylindrical  spreading,  the  second  and 
third  terms  account  for  variations  in  the  medium  density  and  average  angle  of 
propagation  of  the  ray  tube.  The  last  term  shows  the  effects  of  ray  spreading  or 
converging,  which  dominates  the  2D  field  structure. 

The  equation  for  the  evolution  of  the  ray  amplitude,  Equation  (4.21),  reveals  two 
important  features  in  ray  propagation: 


1)  Sz— >0:  when  the  vertical  spread  of  the  rays  Sz  and,  consequently,  the 

area  S  tends  to  zero,  a  high  intensity  region  is  formed  where  all  the  rays  in 
the  tube  converge  to  a  focal  point  or  caustic.  In  this  region,  the  amplitude 
tends  to  infinity,  which  is  a  limitation  of  classical  ray  theory; 


2)  8z<  0:  after  passing  through  a  caustic,  8z  changes  sign  and  the 

amplitude  can  become  complex;  therefore,  Equation  (4.21)  may  be 
rewritten  as 


A(r)  =  A(ro) 


jr„  pc  8(f>n  sin  0o  r0S90 
r  p0c0  8(f)  sin  9  8z 


jN - 


(4.24) 


The  exponential  term  in  Equation  (4.24)  indicates  a  n/l  phase  change  at  each  caustic 
(N  represents  the  number  of  caustics).  Therefore,  to  obtain  the  correct  phase  of  the 
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pressure  field,  it  is  necessary  to  keep  track  of  the  number  of  caustics  encountered 
by  the  rays. 

In  addition  to  the  phase  changes  that  occur  at  caustics,  rays  also  undergo  phase 
changes  at  reflections.  For  reflections  from  the  surface  (which  can  be  treated  as  a 
pressure  release  boundary),  each  ray  picks  up  a  n  phase  shift.  Reflections  at  the  seafloor 
can  become  more  complicated  because  the  phase  change  depends  upon  the  reflection 
angle,  which  can  vary  considerably  among  all  the  rays  contributing  to  the  field  solution  at 
a  point.  Bottom  reflections  are  generally  not  treated  as  ideal,  but  contribute  to  some  loss 
in  amplitude  along  the  ray  path.  The  complete  description  of  the  field  at  a  point  using  a 
ray  approach  may  then  be  defined  as  [30] 

M 

*)=ZMr)l 

m= 1 

where  m  =  1, . . .  M  are  the  number  of  eigenrays,  \Am\  is  the  amplitude  of  the  mth  eigenray, 

9? mn(r'mn}  is  the  reflection  coefficient  (which  may  be  complex)  for  circumstances  in 

which  the  mth  eigenray  encountered  Nbm  bottom  reflections  at  points  r'  ,  Ncm  is  the 

number  of  caustics,  and  Ns  m  is  the  number  of  surface  reflections  of  the  mth  eigenray.  The 

actual  trajectory  of  the  eigenray  is  defined  by  the  path  sm.  Keeping  track  of  all  of  the 
various  factors  that  determine  the  pressure  field  from  ray  calculations  can  be  a 
challenging  task,  and  is  at  the  root  of  all  acoustic  ray  tracing  models. 

e.  Geometric  and  Gaussian  Beams 

In  addition  to  caustics,  the  standard  ray  tracing  method  produces  perfect  shadow 
zones,  which  represent  limitations  of  the  classical  theory.  One  way  to  deal  with  such 
artifacts  is  the  use  of  geometric  or  Gaussian  beams  [26],  [32],  and  [33].  This  approach 
permits  some  leakage  energy  in  the  shadow  zones  and  smooths  out  the  caustics, 
providing  a  good  agreement  with  more  computationally  intensive  full-wave  models. 

As  described  in  [26],  a  beam  is  constructed  around  each  ray  where  the  amplitude 
of  the  beam  is  tapered  to  vary  from  the  central  ray,  decaying  to  zero  on  either  side.  The 


I  I , ) 

n= 1 


-id)  [  m  — — -ds'm 
,  Jo  e(sm) 


elNc,m  2  elNs„ 


(4.25) 
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half-width  W (.v)  of  the  beam  is  defined  by  the  distance  from  the  central  ray  where  the 

beam  amplitude  is  made  to  vanish.  This  width  is  chosen  precisely  so  that  the  beam 
vanishes  at  the  location  of  its  neighboring  ray.  Thus,  the  pressure  for  the  beam  is  given  by 

Pbeam  (s,  n )  =  Abeam  (5)  <j){s,  n )  eit0T{s) ,  (4.26) 

where  (f>{s)  can  be  the  hat-shaped  function  (geometric  beam)  [32]  or  a  Gaussian  function 
(Gaussian  beam)  [34],  respectively 


or 


</>(s,n) 


' W(s)-n 

<  W(s) 


for  n<W (5) 


0  elsewhere 

V 


(4.27) 


(4.28) 


where  n  =  0  represents  the  central  ray.  In  this  approach,  the  ray  becomes  a  central  ray  of 
a  Gaussian  varying  beam,  Equation  (4.28),  or  a  linear  varying  beam,  Equation  (4.27),  as 
depicted  in  Figure  17. 


-IV  (.S')  W(s)  " 


Figure  17.  Geometric  Beams  around  Each  Ray.  Adapted  from  [32], 
f.  BELLHOP  Beam  Tracing  Model 

In  this  dissertation,  the  ray  tracing  code  used  is  BELLHOP,  developed  by  Michael 

Porter  and  available  under  GNU  public  license  in  the  Ocean  Acoustic  Library  website. 

According  to  [24],  BELLHOP  is  a  beam  tracing  model  for  predicting  acoustic  pressure 
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fields  in  ocean  environments  where  several  types  of  beams  may  be  implemented, 
including  Gaussian  and  hat-shaped  beams,  with  both  geometric  and  physics-based 
spreading  laws.  BELLHOP  can  generate  a  variety  of  outputs  including  transmission  loss 
and  eigenrays.  It  allows  for  range-dependence  in  the  top  and  bottom  boundaries 
(altimetry  and  bathymetry),  as  well  as  in  the  sound  speed  profile.  Additional  input  files 
permit  the  specification  of  directional  sources  as  well  as  geoacoustic  properties  for  the 
bounding  media.  Top  and  bottom  reflection  coefficients  may  also  be  included. 

The  eigenrays  (amplitude  and  phase)  represent  the  ideal  channel  impulse  response 
(IR).  This  computation  is  done  exactly  the  same  way  as  is  done  for  a  regular  ray  trace. 
However,  BELLHOP  only  saves  the  rays  whose  associated  beams  make  a  contribution  to 
the  specified  receiver  location. 

2.  Matched  Peak  Impulse  Response  Processor 

In  this  section  an  algorithm  that  makes  use  of  a  ray  tracing  code  to  model  the 
environment  is  developed  to  estimate  the  range  between  a  UUV  and  a  reference  point. 
The  use  of  a  ray  tracing  code  permits  us  to  take  into  account  the  multipath  and  the 
refraction  of  the  sound  waves  in  the  range  estimation. 

a.  Motivation 

Problems  involving  estimates  of  distances  underwater  belong  to  a  broader  class  of 
source  localization  problems  which  have  been  extensively  studied.  The  state-of-the-art 
approach  makes  use  of  techniques  such  as  matched  field  processing  [35]  or  matched  IR 
processing  [36]  to  estimate  several  acoustic  parameters,  including  source  position  and 
horizontal  distances  from  source  to  target. 

Those  techniques  involve  acoustic  field  measurements  (in  the  case  of  the  matched 
field  processing)  or  channel  IR  measurements  (in  the  case  of  the  matched  IR  processing), 
and  acoustic  propagation  models  to  generate  synthetic  patterns.  These  synthetic  patterns 
are  then,  iteratively,  used  to  match  the  complete  measured  pattern. 

This  approach  can  produce  good  results  in  estimating  distances  underwater  but 
demands  very  high  computational  efforts  and  may  take  hours  to  converge.  Therefore,  it  is 
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not  suitable  for  use  in  real-time  applications.  In  what  follows,  it  will  be  shown  that  the 
acoustic  modem’s  matched  filter  output  may  be  considered  for  estimation  of  the  channel 
IR  (measured  estimation).  This  estimation  can  be  emulated,  generating  a  prediction  for 
the  channel  IR  (synthetic  prediction). 

The  synthetic  predicted  channel  IR  can  then  be  used  iteratively  to  match  the 
measured  estimated  channel  IR.  However,  instead  of  trying  to  match  the  entire  time 
series,  as  in  the  traditional  approach,  only  the  arrival  with  the  highest  amplitude  is 
utilized  (see  Section  A).  This  algorithm  shall  run  inside  the  UUV  and  is  fast  enough  to 
work  in  real-time. 

b.  Acoustic  Channel  IR  Estimation 

In  practical  applications  in  which  a  single  source  transmits  to  a  receiver,  an 
estimate  of  the  acoustic  channel  IR  is  obtained  by  transmitting  a  broadband  pulse  (like  a 
chirp)  from  the  source,  which  is  then  match  filtered  (or  cross-correlated)  at  the  receiver. 
Assuming  that  the  received  signal  x'(t )  is  cross-correlated  with  a  replica  of  the 

transmitted  signal,  x(t ),  such  that 

+00 

h(t)=  |  x'{r)x*  {t  +  r^dr ,  (4.29) 

—00 

then  h[t)  may  represent  an  estimation  for  the  channel  IR. 

The  received  signal  can  be  understood  as  being  the  result  of  the  convolution 
between  the  transmitted  signal  and  the  acoustic  channel  complex  IR  plus  noise  (for 
simplicity,  no  attenuation  or  other  propagation  effects  are  being  considered).  Specifically, 

+CX> 

v'(r)  =  1  h[a)x(r-a)da  +  n(r),  (4.30) 

—00 

where  h(a) is  the  true  acoustic  channel  IR.  Substituting  Equation  (4.30)  into  (4.29) 
yields 


43 


-|-UU  -|-UU 

h(t)=  |  J  /z(«)  x(t-cc) da  +  n(j) 


—00  00 
+00  I  +00 


X 


(t  +  r)  ch 


(4.31) 


=  JM  a)\  |  x(t-cc)x*  (t  +  r) dr  \da+  J  n(r)x*  (r  +  r) dr . 


If  the  transmitted  signal  autocorrelation  has  an  impulse-like  behavior,  the  term  in 
brackets  in  Equation  (4.31)  can  be  approximated  as  [36] 

+GO 

^  x{t -a)x*  (t  +  r)dr  &  Ex8(t-a) ,  (4.32) 

—oo 

where  Ex  is  the  energy  of  the  signal. 

Normalizing  Equation  (4.32)  to  the  energy  of  the  signal,  substituting  into 
Equation  (4.31),  and  considering  that  signal  and  noise  are  uncorrelated,  produces 

+oo 

h(t)  ~  |  /*(a)  (J(t-a) da  =  h(t) .  (4.33) 


From  Equations  (4.32)  and  (4.33),  we  may  conclude  that  when  the  transmitted  signal 
autocorrelation  approaches  an  impulse-like  response,  then  the  matched  filter  (or  cross¬ 
correlation)  output  provides  an  accurate  estimate  of  the  channel  IR. 


As  described  in  Chapter  III,  Section  C,  one  of  the  most  common  waveforms  used 
in  underwater  acoustics  is  the  HFM  pulse.  The  half-width  of  the  main  lobe  in  its 
autocorrelation  is  approximately  \IB  (B  is  the  swept  bandwidth),  which  means  that 
increasing  the  transmitted  bandwidth  improves  the  quality  of  the  estimation  for  the 
channel  IR.  The  drawback  of  increasing  the  bandwidth  is  that,  in  order  to  maintain  the 
same  SNR  at  the  receiver,  the  power  of  the  signal  must  be  increased  since  noise  power 
increases  with  bandwidth. 


c.  Matched  Peak  IR  Processor  Algorithm 

Considering  a  single  source  and  receiver,  the  formulation  of  the  ray  theory  in  the 
time  domain  permits  a  fast  assessment  of  the  ideal  acoustic  channel  complex  IR 
(amplitude  and  phase  of  the  eigenrays),  where  each  eigenray  may  be  represented  by  a 
complex  amplitude  function  as 

A(/)  =  j\em ,  i  =  \,2 , ...  M  ,  (4.34) 
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where  A  is  the  amplitude  and  (pt  is  the  phase  of  the  ith  eigenray,  and  M  is  the  number  of 
eigenrays.  Note  that  phase  <pi  accounts  for  boundary  reflections  and  caustics. 

A  synthetic  received  signal  may  be  approximated  by  phased  addition  of  the 
complex  amplitude  associated  with  each  eigenray,  which  resembles  the  convolution 
between  the  ideal  channel  IR  and  a  replica  of  the  transmitted  signal.  According  to  [26], 
Chapter  8,  a  synthetic  received  signal  may  be  constructed  as 

M 

X' '  (0  =  Z \  • X  (*  - Tm  ) -  ■ AlJ  {t - Cn  )  >  (4-35) 

m= 1 

where  \  ,  Ar  ,  and  rm  are,  respectively,  the  real  component  of  the  complex  amplitude, 
the  imaginary  component  of  the  complex  amplitude,  and  travel  time  of  the  mth  eigenray. 
M  is  the  number  of  contributing  eigenrays,  x(t)  is  the  transmitted  signal,  and  x(t)  is  the 
Hilbert  transform  of  the  transmitted  signal. 

To  emulate  the  signal  processing  made  by  the  acoustic  modems,  the  synthetic 
received  signal  must  be  shifted  to  baseband  (for  more  efficient  processing)  and  match 
filtered  considering  a  baseband  replica  of  the  transmitted  signal.  The  output  of  the 
matched  filter  will  represent  the  synthetic  channel  IR.  In  this  work,  the  time  of  the 
highest  peak  is  defined  as  the  predicted  time,  tp. 

In  Section  A,  we  have  seen  that  the  measured  one-way  travel  time,  tm ,  calculated 

by  the  Benthos  acoustic  modems,  is  based  on  the  time  of  the  highest  peak  in  the  matched 
filter  output.  The  approach  used  here  is  to  iteratively  attempt  to  match  (within  a 
tolerance)  tp  with  tm  by  varying  the  horizontal  distance,  R,  between  the  UUV  and  the 
reference  point  in  the  ray  tracing  code.  The  best  estimate  for  R  is  achieved  when 
tm—s  <tp<tm+s ,  where  s  is  the  assumed  tolerance. 

This  approach  is  fast  enough  to  run  on  dedicated  hardware  in  real  time  and  is 
more  accurate  for  estimation  of  range  than  a  simple  multiplication  of  the  one-way  travel 
time  with  the  characteristic  medium  sound  speed.  Any  eventual  errors  are  small  enough 
to  be  handled  by  the  tracking  algorithm,  as  is  described  in  Chapter  V. 
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The  complete  algorithm  is  depicted  in  Figure  18.  Consider  a  UUV  equipped  with 
an  acoustic  modem  operating  in  a  known  region  (bottom  properties  and  bathymetry  are 
known),  which  is  able  to  measure  the  sound  speed  of  the  water  column.  In  this  area, 
assume  that  a  network  of  reference  points  is  available,  as  described  in  Chapter  II. 

The  UUV  and  one  of  the  reference  points  establish  acoustic  communication,  and 
the  acoustic  wave  travel  time  tm  (one-way  travel  time)  is  measured.  As  part  of  this 
measurement,  the  UUV  receives  the  coordinates  of  the  reference  point  (latitude, 
longitude,  and  depth).  Latitude  and  longitude  are  later  used  in  the  tracking  algorithm 
(Chapter  V). 

After  a  successful  travel  time  measurement,  the  algorithm  (running  inside  the 
UUV)  starts  loading  into  the  ray  tracing  code  (BELLHOP)  the  first  guess  for  the 
horizontal  distance  between  the  UUV  and  the  reference  point  (R  =  1500 x  tm ),  having 

already  pre-loaded  the  bottom  properties,  the  bathymetry,  and  measured  sound  speed 
profile.  Then  the  depth  of  the  UUV  and  the  depth  of  the  reference  point  are  loaded  into 
BELLHOP. 

BELLHOP  is  set  to  give  as  output  the  complex  IR  (eigenray  amplitudes  and 
phases)  of  the  ocean  waveguide.  Now,  the  synthetic  received  signal  can  be  constructed 
according  to  Equation  (4.35),  shifted  to  baseband,  and  match  filtered.  The  time  of  the 
highest  peak,  t  ,  is  taken  and  compared  with  the  measured  one-way  travel  time,  tm .  This 

process  is  repeated,  adjusting  the  horizontal  distance,  R,  until  t  and  tm  match  within  the 
tolerance  set  by  s  . 
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Figure  18.  Matched  Peak  IR  Processing  Flowchart 
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d.  Examples  of  Application 

During  the  sea  trials  that  took  place  in  Monterey  Bay  on  May  20,  2015,  several 
successful  travel  time  measurements  between  the  UUV  and  the  command  ship  were 
taken.  In  addition,  the  probe  signal’s  matched  filter  output  associated  with  each 
measurement  was  also  recorded. 

The  area  of  operation,  the  bottom  properties,  the  sound  speed  profile,  and  the 
assumption  made  for  the  beam  patterns  are  described  as  follows. 

The  assets  were  operating  in  a  relatively  flat  area,  around  75  m  of  depth,  at  the 
edge  of  the  Monterey  Bay  submarine  canyon  (see  Figure  19). 
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Figure  19.  Area  of  Operation  for  May  20,  2015,  and  August  12,  2015,  Sea  Trials 


The  bottom  was  modeled  following  Hamilton’s  approach  [37],  and  the  sediments 
were  identified  as  being  silt  clay/sand-silt-clay  with  the  following  characteristics: 

Compressional  sound  speed:  1560  m/s; 

Density:  1.6  g/cm3; 

Attenuation  for  compressional  waves:  5  dB/m. 
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The  sound  speed  profile  was  measured  by  the  UUV  and  its  average  is  depicted  in  Figure 

20. 


Figure  20.  Sound  Speed  Profile  Measured  by  the  UUV  on  May  20,  2015 


For  the  purpose  of  calculations  in  this  dissertation,  the  acoustic  modem’s 
horizontal  and  vertical  transducer  beam  patterns  are  considered  omnidirectional.  This  is 
known  to  be  a  simplification,  and  further  efforts  to  characterize  the  beam  patterns  of  the 
transducers  mounted  in  the  assets  are  planned. 

Example  1:  Mission-4:  20:55:10  (hh:mm:ss)  GMT 

-  UUV  depth:  44  m 

Command  ship  transducer  depth:  5  m 
Measured  one-way  travel  time  (tm):  528.13  ms 
Slant  range  rough  estimate  (1500  x  tm):  792.2  m 
Horizontal  range  rough  estimate:  791.2  m 

Figure  21  shows  (a)  the  paths  of  the  eigenrays  and  (b)  the  amplitude  and  time  of 
arrival  associated  with  these  rays.  Note  that  the  first  group  of  arrivals  (representing  the 
direct  path  and  a  surface  bounce)  has  smaller  amplitude  than  the  second  group  of  arrivals 
(representing  a  bottom  bounce  and  surface-and-bottom  bounce).  This  counterintuitive 
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behavior  is  due  to  the  spreading  of  the  ray  tubes  along  their  trajectories,  according  to  the 
following  explanation. 


Eigenray  before  Matched  IR  peak 


x  ip"3  Ideal  IR  before  Matched  IR  peak,  Sd=5m,  Rd=44m,  R=791.2 


Figure  21.  (a)  Eigenrays  and  (b)  Ideal  IR  before  Matched  Peak  IR  Processing. 

Data  from  Example  1 

There  are  two  rays  for  each  arrival.  Each  one  represents  the  central  ray  of  the 
limiting  tube  that  is  contributing  to  the  pressure  field  at  the  receiver  (Sections  B.l.e  and 
B.l.f).  These  ray  tubes  undergo  spreading  and  shrinking  along  their  trajectories.  The  rate 
of  spreading  or  shrinking  is  maximized  in  regions  of  large  sound  speed  gradients.  Note 
that  in  the  absence  of  a  sound  speed  gradient,  the  rays  travel  in  straight  lines  and  the  rate 
of  ray  spreading  is  constant  with  range.  For  the  sound  speed  in  this  environment,  the 
sound  speed  gradient  is  primarily  negative  at  all  depths  (see  Figure  20),  but  the  largest 
gradients  are  near  the  surface  where  the  ocean  temperature  varies  significantly  with 
depth.  Let  us  consider  two  situations  (referring  to  Figure  21): 

-  Those  rays  that  arrive  first  follow  a  direct  path  and  a  nearby  path  that 
reflects  from  the  surface.  The  entire  propagation  path  of  these  rays  is 
concentrated  in  the  shallowest  part  of  the  water  column  where  the  sound 
speed  gradient  is  maximum.  Thus,  the  rate  at  which  the  associated  ray 
tubes  spread  is  also  maximized,  thereby  reducing  the  amplitude  observed 
at  the  receiver. 
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-  The  next  set  of  ray  arrivals  corresponds  to  steeper  propagation  paths  that 
do  not  remain  in  the  upper  part  of  the  water  column.  Instead,  most  of  their 
propagation  path  is  concentrated  in  the  lower  regions  of  the  water  column 
where  the  sound  speed  gradient  is  significantly  smaller.  The  result  is  that 
the  cumulative  ray  tube  spreading  along  these  paths  is  less  than  the  direct 
path,  and  consequently  the  contribution  to  the  pressure  field  at  the  receiver 
is  larger,  despite  the  longer  path  lengths  and  losses  at  the  boundaries. 

Figure  22  provides  a  comparison  between  (a)  the  acoustic  modem  matched  filter 
output  (measured  channel  IR  estimation — real  data)  and  (b)  a  matched  filter  output 
prediction,  as  described  in  Section  B.2.c  (predicted  channel  IR).  To  save  memory,  the 
modem  did  not  record  the  complete  matched  filter  output  time  series.  For  this  reason,  the 
time  axes  are  different  on  the  plots.  Therefore,  the  time  of  the  highest  peak  in  Figures 
22a,  25a,  and  28a  should  be  understood  as  being  the  measured  one-way  travel  time  (tm) 
instead  of  what  is  depicted.  Nevertheless,  the  numbers  indicated  on  Figures  22a,  25a,  and 
28a  are  still  valid  to  calculate  the  time  difference  between  the  highest  peak  and  the  first 
arrival. 
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Acoustic  Modem  Matched  Filter  Output,  Mission  4  May  20,  2015  20:55:10  GMT 


Figure  22.  Acoustic  Modem  (a)  Matched  Filter  Output  and 
(b)  Matched  Peak  IR  Algorithm  Output. 
Data  from  Example  1 


In  Section  B.2.c,  an  algorithm  to  match  the  time  of  the  highest  peak  in  the 

predicted  matched  filter  output  (predicted  channel  IR)  with  the  measured  one-way  travel 

time  was  developed  making  use  of  a  ray  tracing  code  to  model  the  propagation  through 

the  environment.  Figure  22b  represents  the  result  of  this  algorithm  after  minimizing  the 
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difference  between  measured  and  predicted  highest  peak  times.  Note  that  the  time  of  the 
highest  peak  in  Figure  22b  matches  with  the  one-way  travel  time  measured  by  the 
acoustic  modem. 

Another  important  point  is  that  there  is  a  match  between  real  data  (Figure  22a) 
and  the  prediction  (Figure  22b)  when  the  time  difference  between  the  first  and  the  highest 
peaks  is  compared.  Probably  due  to  noise  or  some  kind  of  scattering  (in  the  water  column 
or  at  the  boundaries),  no  good  match  in  the  last  portion  of  Figure  22  (a)  and  Figure  22  (b) 
can  be  found. 

The  final  result  depicted  in  Figure  23  shows  that  the  horizontal  distance  estimated 
when  the  peaks  are  matched  is  782.4  m,  compared  with  the  original  rough  estimate  of 
791.2  m.  The  range  correction  provided  by  the  algorithm  is  8.8  m,  in  this  case. 


Eigenray  after  Matched  IR  peak 


x  10"3  Ideal  IR  after  Matched  IR  peak,  Sd=5m,  Rd=44m,  R=782.38 


Figure  23.  (a)  Eigenrays  and  (b)  Ideal  IR  after  Matched  Peak  IR  Processing. 

Data  from  Example  1 


Example  2:  Mission-2:  19:04:41  (hh:mm:ss)  GMT 

-  UUV  depth:  1.5  m 

Command  ship  transducer  depth:  5  m 
Measured  one-way  travel  time  (tm):  323.6  ms 
Slant  range  rough  estimate  (1500  x  tm):  485.4  m 
Horizontal  range  rough  estimate:  485.4  m 

Figure  24  shows  (a)  the  eigenray  paths  and  (b)  the  amplitude  and  time  of  arrival 


associated  with  these  rays.  Note  that  the  first  group  of  arrivals  (representing  the  direct 
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path  and  a  surface  bounce)  has  greater  amplitude  when  compared  with  the  second  group 
of  arrivals  (representing  a  surface-and-bottom  bounce  and  surface-bottom-and-surface 
bounce). 


Figure  24.  (a)  Eigenrays  and  (b)  Ideal  IR  before  Matched  Peak  IR  Processing. 

Data  from  Example  2 


In  this  case,  the  attenuation  and  bottom  losses  associated  with  the  longer  path 
length  are  the  major  factors  that  cause  the  smaller  amplitude  in  the  second  group  of 
arrivals  (Figure  24b).  Another  interesting  point  is  the  small  time  (and  path  length) 
difference  between  the  direct  path  and  the  surface  bounce  in  the  first  group  of  arrivals. 
Due  to  this  small  time  difference,  surface  interference  may  play  an  important  role  in  the 
result  presented  by  the  matched  filter  output. 

As  before,  the  time  of  the  highest  peak  in  Figure  25b  matches  with  the  one-way 
travel  time  measured  by  the  acoustic  modem,  and  there  is  a  match  between  real  data 
(Figure  25  a)  and  the  prediction  (Figure  25b)  when  the  time  difference  between  the  first 
and  the  highest  peaks  is  compared.  Again,  probably  due  to  noise  and  some  kind  of 
scattering,  no  match  in  the  last  portion  of  Figure  25  (a)  and  Figure  25  (b)  can  be  found. 
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Figure  25.  Acoustic  Modem  (a)  Matched  Filter  Output  and 
(b)  Matched  Peak  IR  Algorithm  Output. 
Data  from  Example  2 


It  is  interesting  to  note  the  effect  of  surface  interference  in  the  first  group  of 
arrivals.  Due  to  the  small  time  difference,  the  surfaced  reflected  wave  (180°  out  of  phase) 
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almost  completely  cancels  the  direct  wave  contribution,  causing  a  reduction  in  the 
matched  filter  output  for  the  first  peak  (see  Figure  25b). 

The  measured  data  shows  similar  behavior  (Figure  25a)  although,  apparently,  the 
reduction  in  the  first  peak  was  not  as  severe  as  predicted  by  the  ray  model.  An  attenuation 
in  the  second  group  of  arrivals  associated  with  the  UUV  and  command  ship  transducer’s 
beam  patterns  may  be  one  contributing  factor  for  this  behavior;  the  beam  patterns  were 
considered  omnidirectional  in  all  calculations.  Changes  in  the  transducer’s  depth,  due  to 
wave  motion,  in  both  assets  may  also  be  a  contributing  factor,  as  could  scattering  of  the 
surface  reflected  path.  In  this  particular  case,  the  UUV  was  at  the  surface  in  the  “comms” 
position  (i.e.,  with  head  tilted  down  and  tail  [the  antenna]  tilted  up),  which  can  be  heavily 
affected  by  surface  wave  motion.  The  changing  in  depth,  caused  by  the  wave  motion, 
may  affect  the  amplitude  of  the  arrivals. 

The  final  result,  depicted  in  Figure  26,  shows  that  the  horizontal  distance 
estimated  when  the  peaks  are  matched  is  461.7  m,  compared  with  the  original  rough 
estimate  of  485.4  m.  In  this  case,  the  range  correction  provided  by  the  algorithm 
is  23.7  m. 


Eigenray  after  Matched  IR  peak 


Range  (m) 


Figure  26.  (a)  Eigenrays  and  (b)  Ideal  IR  after  Matched  Peak  IR  Processing. 

Data  from  Example  2 
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Example  3:  Mission-4:  20:41:06  (hh:mm:ss)  GMT. 

-  UUV  depth:  48.2  m 

Command  ship  transducer  depth:  5  m 
Measured  one-way  travel  time  (tm):  430  ms 
Slant  range  rough  estimate  (1500  x  tm):  645  m 
Horizontal  range  rough  estimate:  643.6  m 

Figure  27  shows  (a)  the  eigenray  paths  and  (b)  the  amplitude  and  time  of  arrival 
associated  with  these  rays.  Note  that  the  first  group  of  arrivals  (representing  the  direct 
path  and  a  surface  bounce)  has  smaller  amplitude  when  compared  with  the  second  group 
of  arrivals  (representing  a  bottom  bounce  and  surface-and-bottom  bounce).  Again,  this  is 
due  to  the  spreading  of  the  ray  tubes  along  their  trajectories. 


Eigenray  before  Matched  IR  peak 


x  io"3  Ideal  IR  before  Matched  IR  peak,  Sd=5m,  Rd=48.2m,  R=643.6 


Time  (s) 


Figure  27.  (a)  Eigenrays  and  (b)  Ideal  IR  before  Matched  Peak  IR  Processing. 

Data  from  Example  3 


As  in  the  previous  examples,  the  time  of  the  highest  peak  in  Figure  28b  matches 
with  the  one-way  travel  time  measured  by  the  acoustic  modem,  and  there  is  a  match 
between  real  data  (Figure  28a)  and  the  prediction  (Figure  28b)  when  the  time  difference 
between  the  first  and  second  peaks  is  compared.  Again,  probably  due  to  noise  and  some 
kind  of  scattering,  it  is  hard  to  find  a  good  match  in  the  last  portion  of  Figure  28  (a)  and 
Figure  28  (b). 
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Acoustic  Modem  Matched  Filter  Output,  Mission  4  May  20,  2015  20:41:06  GMT 


Matched  IR  Peak  Algorithm  Output 


Figure  28.  Acoustic  Modem  (a)  Matched  Filter  Output  and 
(b)  Matched  Peak  IR  Algorithm  Output. 
Data  from  Example  3 


In  Figure  27,  the  second  group  of  arrivals  (representing  a  bottom  bounce  and 
surface-and-bottom  bounce)  has  higher  amplitude  than  the  first  group  of  arrivals 
(representing  the  direct  path  and  a  surface  bounce).  However,  Figure  28b  indicates  that 
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the  interference  between  arrivals  produced  the  opposite  effect,  and  the  first  group 
of  arrivals  produces  the  highest  IR  peak.  This  effect  was  confirmed  in  the  real  data 
(Figure  28a). 

Also  observable  in  the  Figure  28  plots  is  a  mismatch  of  the  amplitude  difference 
between  the  first  and  second  arrival  groups  in  the  real  data  compared  with  the  prediction. 
This  mismatch  could  be  associated  with  the  transducer  beam  patterns  on  the  UUV  and 
command  ship,  surface  scattering,  or  unknown  environmental  effects. 

The  final  result,  depicted  in  Figure  29,  shows  that  the  horizontal  distance 
estimated  when  the  peaks  are  matched  is  642.2  m,  compared  with  the  original  rough 
estimate  of  643.6  m.  The  range  correction  provided  by  the  algorithm  is  only  1.4  m.  This 
small  difference  simply  accounts  for  the  refraction  of  the  sound  waves  (rays  bending 
versus  straight  line  propagation),  since  it  is  this  first  arrival  that  provides  the  highest  peak 
in  the  matched  filter  output. 


Eigenray  after  Matched  IR  peak 


x  io'3  ldeal  IR  after  Matched  IR  peak,  Sd=5m,  Rd=48.2m,  R=642.18 


Time  (s) 


Figure  29.  (a)  Eigenrays  and  (b)  Ideal  IR  after  Matched  Peak  IR  Processing. 

Data  from  Example  3 

To  claim  any  improvement  associated  with  the  method  developed  in  this  chapter, 
we  use  the  horizontal  distances  estimated  when  the  peaks  are  matched  and  roughly 
estimated  as  inputs  for  the  tracking  algorithm,  as  developed  in  Chapter  V,  and  the  results 
are  compared  in  Chapter  VII. 
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V.  THE  TRACKING  ALGORITHM 


Parts  of  this  chapter  were  previously  published  by  EEEE  [1].* 

As  defined  in  [38],  “estimation  is  the  process  of  inferring  the  value  of  a  quantity 
of  interest  from  indirect,  inaccurate  and  uncertain  observations,”  and  “tracking  is  the 
estimation  of  the  state  of  a  moving  object  that  can  be  done  using  one  or  more  sensors  at 
known  locations.”  The  Kalman  filter,  in  all  its  variations,  is  one  of  the  most  widely  used 
methods  for  tracking  and  estimation. 

In  this  chapter,  a  tracking  algorithm  based  on  the  KF  is  developed  for  a  UUV. 
Due  to  non-linearities  in  the  measurement  equation  two  approaches  are  taken,  first  using 
the  EKF  and  then  using  the  UKF.  The  advantages  and  disadvantages  of  both  approaches 
are  highlighted.  Additionally,  tuning  and  smoothing  algorithms  are  discussed  and,  as  sea 
current  is  part  of  the  UUV  state  representation,  a  consensus  algorithm  to  take  advantage 
of  these  predictions  from  different  UUVs  is  developed. 

A.  KALMAN  FILTER 

The  KF  [39]  can  be  seen  as  an  optimal  state  estimator  for  discrete-time  linear 
stochastic  dynamic  systems  for  which  state  space  representation  is  given  by  Equations 
(2.13)  and  (2.14).  KF  is  the  optimal  minimum  mean  square  error  estimator  when  all  the 
noises  entering  the  system  (including  the  initial  state  error)  are  Gaussian.  If  these  random 
variables  are  not  Gaussian,  KF  is  still  the  best  (minimum  error  variance)  linear  state 
estimator  [38], 


*  ©  2016  IEEE.  Personal  use  of  this  material  is  permitted.  Permission  from  IEEE  must  be  obtained  for 
all  other  uses,  in  any  current  or  future  media,  including  reprinting/republishing  this  material  for  advertising 
or  promotional  purposes,  creating  new  collective  works,  for  resale  or  redistribution  to  servers  or  lists,  or 
reuse  of  any  copyrighted  component  of  this  work  in  other  works. 
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KF  is  a  recursive  filter  composed  of  two  steps: 

-  Prediction:  the  state  of  the  system  is  predicted  according  to  the  state 
equation,  given  the  previous  updated  state  estimate; 

-  Measurement  update:  prediction  is  updated  given  a  measurement  collected 
at  this  time  step. 

Derivations  of  the  KF  equations  can  be  found  in  [9],  [38]  and  are  not  presented  in  this 
dissertation. 


a)  System  model  and  initial  assumptions 


The  state  and  measurement  equations  for  a  discrete-time  linear  stochastic  dynamic 
system,  as  previously  stated  in  Equations  (2.13)  and  (2.14),  are 

x(k  +  l)  =  F(k)x(k)  +  G(k)u(k)  +  o(A:) 

z(k)  =  K(k)x(k)  +  <o(k),  (5-1) 

where,  for  simplicity,  D  in  Equation  (2.14)  was  set  to  zero,  and  F,  G,  u,  and  H  are 
deterministic  sequences  assumed  to  be  known  for  all  k.  The  noises  o(/c)  and  co(k)  are 

assumed  to  be  zero  mean,  white  (i.e.,  uncorrelated)  Gaussian  processes,  such  that 


£[o(k)]  =  £[©(£)]  =  0 

£[o(t)o(y)r]  =  QAj 

£[<o(*r)ra(y) 

E|  v(k)(o(j)T 


=  R  A; 


■OVkJ. 


(5.2) 


At  each  step  ( k ),  the  state  estimate  and  error  covariance  matrix  is  determined  by  a 
sequence  of  well-defined  prediction  and  measurement  updates  (also  called  correction),  as 
follows: 


b)  Prediction 

x(k  +  l\k)  =  ¥(k)x(k\k)+G(k)u(k) 
F(k  +  l\k)  =  F(k)F(k\k)F(k)T  +Qk , 
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where  x(k|k)  is  the  last  updated  state  estimate  and  P(&|&)  is  the  covariance  matrix 
associated  with  this  estimate. 

c)  Measurement  update: 

Kt+1  =  P(/c  +  l|/c)H(/c  +  l)r  +  l)P(/c  +  l|/c)H(/c  +  l)r  +  R*+1 

x(it  +  l|it  +  l)  =  x(it  +  l|it)  +  Kt+1[z(jk  +  l)-H(jk  +  l)x(jk  +  l|it)]  (5.4) 

P (k  + 1|  k  + 1)  =  [I  - K,+1H  (k  + 1)]  P (k  + 1|  k) [I  -  K,+1H (k  + 1)]"  +Kk+lRk+lKkJ  , 

where 

-  K/;+l  is  called  the  Kalman  Gain  matrix; 

-  x(fc  +  l|fc  +  l)  is  the  updated  state  and  P(k  +  l|k  +  l)  is  the  covariance  matrix 
associated  with  the  updated  state; 

-  z(k  +  l)  is  the  measurement.  If  no  measurement  is  available  at  this  step,  just  set 
x(fc  +  l|fc  +  l)  =  x(fc  +  l|fc)  and  P(k+l|k  +  l)  =  P(k+l|k)  ,  ignoring  step  (c). 

B.  EXTENDED  KALMAN  FILTER 

In  many  cases,  dynamic  systems  are  not  linear  by  nature,  which  means  that  KF 
cannot  be  used  for  state  estimation.  In  such  situations,  the  state  or  the  measurement 
equations  can  be  non-linear.  The  most  common  approach  makes  use  of  Taylor  series 
expansions  to  treat  the  non-linearities,  giving  origin  to  the  EKF. 

a)  System  model  and  initial  assumptions: 

The  state  and  measurement  equations  for  a  discrete-time  non-linear  stochastic 
dynamic  system,  as  previously  stated  in  Equations  (2.11)  and  (2.12),  are 

x  (k  + 1)  =  f  f  x(k)]  +  o  (k) 

r  n  (5-5) 

z(&)  =  h[x(k)]  +  (o(k) , 

where,  for  simplicity,  g[u(k)]  in  Equation  (2.11)  was  set  to  zero  assuming  no 
deterministic  inputs,  and  f[x(k)]  and  h[x(fc)]  are  non-linear  functions.  The  noises 
o(k)and  to(k)  are  assumed  to  be  zero  mean,  white  (i.e.,  uncorrelated)  Gaussian 
processes,  as  in  Equation  (5.2). 
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b)  Taylor  series  expansion: 


In  the  classical  approach,  the  non-linear  functions  are  expanded  to  the  first  order 
around  the  latest  estimate  (second  order  expansions  can  be  found  in  [38]),  as 


f  [x(k)]  =  f  [x(k\k)]  +  Ft  [x(k) - x(k\k)] 
h  [x(k)]  =  h  [x(k  + 1|  k)]  +  H,+1  [x(k)  -  x(k  + 1|  k)] , 


where 


df  [x(fc>] 
dx(k ) 


lx(fc)=x(fc|fc) 


and  H,+1  = 


dh[x(&)] 
d  x(k) 


lx(fc)=x(fc-i-l|&) 


Substituting  Equation  (5.6)  into  Equation  (5.5),  we  obtain 

x(&  +  l)  =  f  [^(^l^iJ  +  Fj,  [x(&)-x(&|&)]  +  o(&) 

=  Fkx(k)+f  [x(&|&)]-Fj,x(&|&)  +  o(&) 
z  (£)  =  h  [x(fc  + 1|  *)]  +  H(t+1  [: x(k )  -  x(k  + 1|  k)]  +  to  (k) 

=  H  k+lx(k) + h  [x(k  + 1|  jfc)]  -  H  Mx(k  + 1|  k)  +  a>(k ) . 


Hence,  we  have  the  linearized  state  and  measurement  equations 

x  (&  + 1)  =  F,,x(£)  +  v  (&)  +  u  (k  ) 
z(t)  =  H 

k+ 1  x{k)  +  y(k)  +  (o(k), 


with  deterministic  terms 


y(ife)  =  f[x(ik|ife)]-Fiti(Jk|Jk) 

y(k)  =  h  [x(k  + 1|  k)~\  -  Hi+1x(£  + 1|  k) . 

c)  Prediction: 

x(k  + 1|  =  F^x^l  k)  +  y(k) 

=  Fj.x(A:|  k) +f  [x(£|  &)]  -  F(.x(^|  k) 
=  {[x(k\k)~\ 


(5.7) 

(5.8) 


(5.9) 


(5.10) 


(5.11) 


P(k  +  l\k)  =  FkP(k\k)FkT  +Qk .  (5.12) 

d)  Measurement  update: 

K,+1=P(/:  +  l|/:)H,+1r[H,+1P(^  +  l|^)H,+1r+R,+1J1 
x(ik  +  l|jk  +  l)  =  x(it  +  l|jfe)  +  Kt+1  [z(ik  +  l)-h[x(it  +  l|ifc)]] 

P  (£  + 1|  £  + 1)  =  [I  -  K(tlHw  ]  P  (fc  + 1|  £ )  [I  -  Kt+1Hi+1  f  +  Kk+lRk+lKk+lT , 
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(5.13) 


where  z(k  + 1)  is  the  measurement.  If  no  measurement  is  available  at  this  step,  just  set 
x(fc+l|fc  +  l)  =  x(fc+l|fc)  and  P(fc  +  l|fc+l)  =  P(fc  +  l|fc),  ignoring  step  (d). 


C.  UNSCENTED  KALMAN  FILTER 


The  UKF  as  seen  in  [40],  [41],  and  [42]  is  a  technique  that  utilizes  the  unscented 
transform  (UT)  [43],  [44],  and  can  be  applied  in  models  of  the  form  of  Equation  (5.5). 
Instead  of  using  linear  approximations,  as  EKF  requires,  the  idea  of  the  UT  is  to 
deterministically  choose  a  fixed  number  of  sigma  points  that  capture  the  mean  (state)  and 
covariance  of  the  original  distribution.  These  sigma  points  are  then  propagated  through 
the  non-linearity,  and  the  mean  and  covariance  of  the  transformed  variable  are  estimated 
from  them.  According  to  [42],  UT  is  able  to  capture  higher  order  moments  when 
compared  with  Taylor  series  based  approximations. 

Assuming  the  system  model  (state  and  measurement  equations)  and  initial 
assumptions  as  in  Equation  (5.5),  at  each  step  ( k ),  the  general  approach  is  along  similar 
lines  as  the  KF  and  the  EKF,  by  prediction  and  measurement  update  as  follows  [45] : 

a)  Prediction 


Form  sigma  points: 


X[0)=x(k\k), 


xi°  =x(k\k^  +  sfn  +  l  <Jp Jk\k) 

xr =mk)-^\sjrW) 


(5.14) 


where  the  matrix  square  root  denotes  a  matrix  such  that  yfPyfp7  =  P ,  n  is  the  length  of  the 
state,  [•].  denotes  the  ith  column  of  the  matrix,  and  A  is  a  scaling  parameter  defined  as 

A  =  a2  (n  +  K)-n  .  The  parameters  a  and  K  determine  the  spread  of  the  sigma  points 
around  the  mean. 


-  Propagate  the  sigma  points  through  the  state  model: 

*  =0,...,2n .  (5.15) 
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-  Compute  the  predicted  state  and  covariance: 

2  n 

i(fc+1l*)  =  Zw/")x*° 

i=0 

2/t  ^ 

p(*+1l*)  =  Zw/c)(x*)-x(*+i|*))(xi0-*(*+i|*))  +Q*, 

1=0 

where  the  weights  W/m)  and  W;.(e)  are  defined  as 


(5.16) 


2 

+  2  ’ 

wT 

-  *  +(.- 
n  +  2 

+/?)> 

W<m) 

1 

i  = 

T  i 

_  2(n  +  2)  ’ 

w<c) 

1 

i  =  l,...,2n, 

l 

2(n  +  2)’ 

(5.17) 


and  /?  is  a  parameter  that  can  be  used  to  incorporate  prior  information  on  the  (non- 
Gaussian)  distribution  of  the  state. 


b)  Measurement  update 
-  Form  the  sigma  points: 

Zw  =  x(fc  +  l|*), 

xi+i =x(fc+i|*)+>/?mr>/p(*+i|fc) 
xL+r =x(*+i|*)->/n+xr>/p(fc+i|*) 


(5.18) 


i  = 


-  Propagate  the  sigma  points  through  the  measurement  model: 

Ci  =  h[*>X*ii]>  i  =  0,...,2n.  (5.19) 

-  Compute  the  predicted  measurement  mean  juk  ,  covariance  of  the  measurement 
Sk ,  and  cross-covariance  of  the  state  and  measurement,  Ck  : 


Hw  -SW,<’>Z® 

1=0 

2  n  j 

St+i  =  Z  (Ci  -  Pw )  (Ci  -  )  +  Rk+ 1 

i=0 

2n  T 

c,+i  =  I  W;<c)  ( -  i  (* + 1|  *))  (C.  -  n*+i )  • 


(5.20) 
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-  Compute  the  Kalman  Gain  and  the  updated  state  and  covariance: 

Ki+1=Ct+IS^ 

i(*  +  l|Jfc  +  l)  =  i(jfc  +  l|ifc)  +  Kt+1  (z(k  +  l)-p,+1)  (5-21> 

P(*  +  l|*  +  l)  =  P(*  +  l|A)-K/t+ISi+1Ki[+,. 

where  z(k  + 1)  is  the  measurement.  If  no  measurement  is  available  at  this  step,  just  set 
x(k  +  l|£;  +  l)  =  x(k  +  l|k)  and  P(k  +  l|k  +  l)  =  P(k  +  l|k),  ignoring  step  (b). 


D.  KF-BASED  TRACKING  ALGORITHMS 

For  convenience  and  completeness,  some  equations  in  this  section  are  repeated 
from  Sections  B  and  C  of  this  chapter,  and  from  Chapter  II,  Section  B.3. 

1.  EKF-Based  UUV  Tracking  Algorithm 

As  demonstrated  in  Chapter  II,  the  UUV  state  space  representation  is  composed 
by  a  linear  state  equation  and  a  non-linear  measurement  equation. 

a)  System  model  and  initial  assumptions: 

From  Equation  (2.28),  the  UUV  state  and  measurement  equations  are 
x(*+l)  =  F(*)x(*)+G(*)i.(*)+o(*) 
z(k)  =  h[k,x(k)]  +  co(k) , 

where 


-  The  matrices  F(k)  and  G(k)  are  deterministic  and  known  for  all  k: 


F  (*) 


"l 

0 

At  (k) 

0 

'A  /(*) 

0 

0 

1 

0 

A  t(k) 

,  G(^)  - 

0 

A  t(k) 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

-  The  input  sequence  u(k)is  a  deterministic  and  known  for  all  k: 


u(k)  =  V(k)cos(^(k)) 


sin(<9(k)) 

cos(#(k)) 


-  z(k)  =  r(k)2  —xr  (k)2  -  yr  (k)2  is  the  measurement; 
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-  h[k,x(k)]  =  x(k)2  -2 x[k)xr  [k)-2y(k)yr  (k)  +  y(k)2  is  a  non-linear  function; 

-  The  noises  o(k) and  co(k)  are  assumed  to  be  zero  mean,  white  (i.e., 
uncorrelated)  Gaussian  processes  as  in  Equation  (5.2). 

In  the  measurement,  r(k )  is  the  horizontal  distance  from  the  UUV  to  the 
reference  point  estimated  according  to  Chapter  IV  (Section  B.2.c).  The  coordinates  of  the 
reference  point  (longitude  and  latitude  mapped  to  Cartesian  coordinates)  are  xr(k )  and 

3VW- 


b)  Taylor  series  expansion 


z  ( k )  =  Hi+1x(k )  +  y  ( k )  +  co  ( k ) 
y  (*)  =  h  [x(k  + 1|  k)]  - H k+1x(k  + 1|  k) , 


(5.23) 


where 


d  h[x(fc)] 
dx(k) 


x(&)=x(k+l|k) 


5h[x(&)]  5h[x(&)]  5h[x(£)]  Sh[x(&)] 

dx(k )  dcx{k)  dcy  (k) 


(5.24) 


x(k)=x(k+l|k) 


=  \2x{k)-2xr{k)  2y(k)-2yr(k)  0  0]^^ 

=  [2x(k  +  l\k)-2xr(k)  2y(k  +  l\k)-2yr(k)  0  o]. 

The  final  result  in  Equation  (5.24)  is  also  known  as  the  Jacobian  matrix, 
c)  Prediction 

*(*+W  =  F(*)*W*)+G(*M*) 

P(k  +  l\k)  =  F(k)P(k\k)F(k)T  +Qk  , 

where  x(k|fc)  is  the  last  updated  state  estimate  and  P(k|fc)  is  the  4x4  covariance  matrix 
associated  with  this  estimate.  Modeling  o(k)  as  continuous  white  noise  acceleration, 
according  to  [38],  the  covariance  matrix  Q,  takes  the  form 


Q(k)  =  cov(u(k))  =  q 2 


Ar(£)3/3  At(k)2 /2  0  0 

Ar(£)3/3  At(kf /2  0  0 

0  0  At(kf  /2  A  t(k) 

0  0  At(kf /2  A  t(k) 


(5.26) 
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where  q  is  the  tuning  parameter  and  At  is  the  sampling  time. 

d)  Measurement  update 

z(k  +  1)  =  r{k  + 1)2  —  xr  [k  +  l)2  —  yr  (k  +  \f 

K,+i=P(^  +  l|/c)H,+1r[H,+1P(/c  +  l|^)H,+1r+R,+1J1  (521) 

x  (jfc  + 1|  k  + 1)  =  x  (k  + 1|  k)  +  Ki+1  [z  (k  + 1)  -h[x(fc  + 1|  *)]] 

P(*  +  l|*  +  l)  =  [l-Ki+1Hi+JP(fc  +  l|fc)[l-K4+1Hi+J%Ki+1Rt+1Ki/, 
where  z  ( /c  +  1)  is  the  measurement.  If  no  measurement  is  available  at  this  step,  just  set 
x(&  +  l|&  +  l)  =  x(&  +  l|&)  and  P(fc  +  l|&  +  l)  =  P(^+l|^),  ignoring  step  (d). 


2.  UKF-Based  UUV  Tracking  Algorithm 

Assuming  the  system  model  and  initial  assumptions  as  in  Equation  (5.22),  at  each 
step  ( k )  the  following  operations  are  performed  [45] : 


a) 


Prediction 


Form  sigma  points: 


i?  =£(*|fc). 


if  =x(k\k)  +  yfc+IyF(k\k) 
Xk+n)  =x(^|^)-V«  +  ^  ^P Jk\k) 


(5.28) 


where  the  matrix  square  root  denotes  a  matrix  such  that  -Jp- \jf*T  =  P ,  n  is  the  length  of 
the  state,  [•].  denotes  the  ith  column  of  the  matrix,  and  X  is  a  scaling  parameter  defined 

as  X  =  a2  (n  +  K)  —  n .  The  parameters  a  and  ic  determine  the  spread  of  the  sigma  points 
around  the  mean.  The  vector  x(A'|  A:  j  is  the  last  updated  state  estimate  and  P(A| A’)  is  the 
4x4  covariance  matrix  associated  with  this  estimate. 

-  Propagate  the  sigma  points  through  the  state  model: 

Xr=F(£)xr+G(£)u(>),  i  =  0,...,2n  ,  (5.29) 

where  F ( A' ) ,  G(A) ,  and  u(A)  are  defined  in  Equation  (5.22). 


69 


-  Compute  the  predicted  state  and  covariance: 


2  n 


'=“  (5.30) 

2  n  j 

p(^  +  1|^)  =  £^)(xf-x(^  +  1|^))(5cf-x(^  +  1|^))  +  Q,, 

(=0 

where  the  matrix  Qk  is  defined  in  Equation  (5.26),  and  the  weights  W,(m)and  W;(e)are 
defined  as 


w<m)  - 

vv0  — 

w<c)  = 

vvo 

W<m>  - 

T  i 

w.(c)  = 


2 

h  -\-  2 


2 

n  +  2 


1 

2(n  +  2)  ’ 

1 

2(n  +  2)  ’ 


i  =  l,...,2n , 

i  =  l,...,2n , 


(5.31) 


and  J3  is  a  parameter  that  can  be  used  to  incorporate  prior  information  on  the  (non- 
Gaussian)  distribution  of  the  state. 


b)  Measurement  update 


Form  the  sigma  points: 

X*+i  =  x(fc  +  l|fc), 

=x(^k  +  l\k^  +  sfn  +  A  ^P(fc  +  l|&) 
:(k  +  \\k)-4n  +  l  ^P(it  +  l|it) 


(5.32) 


x;r=xi 


i  — 


-  Propagate  the  sigma  points  through  the  measurement  model: 

z®=h[l,x^i],  i  =  0,...,2n ,  (5.33) 

where  h^,/^x]  =  x(,)  (&)2  -2x(,)  (&)xr  ( k)-2y(,)  (k)yr  (&)  +  ;y(,)  (£)2 ,  x(,)is  the  position 
in  x  of  the  ith  sigma  point  and  y(,)  is  the  position  in  y  of  the  i,h  sigma  point. 

-  Compute  the  predicted  measurement  mean  /uk  ,  covariance  of  the  measurement 
Sk ,  and  cross-covariance  of  the  state  and  measurement,  Ck : 
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(534) 


“ZW'-’z™ 

i= 0 

2  n  j 

s,„  =ZW,<-’(z<'|1-ft„)(z<"-|i,.1)  +  R„„ 

i=0 

2n  T 

c,+l  =  E  w,(c)  (X«  -  X (*  + 1|  * ))  (z<»  -  pt+1 )  . 

i=0 

-  Compute  the  Kalman  Gain  and  the  updated  state  and  covariance: 

z  (k  + 1)  =  r(k  + 1)2  -  xr  ( k  + 1)2  -  yr  (k  + 1)2 
K,+1=C,+1S^ 

i(*  +  l|t  +  l)  =  x(t  +  l|*)  +  KM(z(*  +  l)-^.,) 

P(*+l|t  +  l)  =  P(t  +  l|*)-K„1S,.,Kl,l , 

where  z( /<  + 1)  is  the  measurement.  If  no  measurement  is  available  at  this  step,  just  set 
x(fc  +  l|fc  +  l)  =  x(fc  +  l|fc)  and  P(^  +  l|^  +  l)  =  P(^  +  l|^),  ignoring  step  (b). 

E.  GENERAL  OBSERVATIONS 

In  this  section,  the  advantages  and  disadvantages  of  using  EKF  and  UKF  are 
highlighted,  and  an  alternative  model  considering  non-additive  Gaussian  noise  is  briefly 
discussed. 


1.  EKF  and  UKF  Comparison 

As  seen  in  Section  B,  EKF  makes  use  of  Taylor  series  expansions  to  treat  the  non- 
linearities  in  the  state  space  model.  This  approach  introduces  limitations  in  the  use  of  this 
technique.  First,  the  non-linearities  need  to  be  differentiable  and  second,  higher  order 
moments  are  not  considered  in  the  Taylor  series  expansion. 

As  presented  in  Section  C,  instead  of  linear  approximations,  UKF  makes  use  of 
deterministically  chosen  sigma  points  to  capture  the  true  mean  and  covariance  of  the 
states.  Those  points  are  then  propagated  through  the  non-linearities,  and  the  mean 
(predicted  state)  and  covariance  are  estimated  from  them.  This  process  is  also  called 
statistical  linearization. 


71 


Advantages  of  the  UKF  in  comparison  with  EKF  are: 

-  Non-linearities  do  not  need  to  be  differentiable  nor  do  their  Jacobian 
matrices  need  to  be  calculated; 

-  UKF  can  capture  higher  order  moments  better  than  EKF  [42]. 

Despite  these  advantages,  UKF  requires  more  computing  efforts  due  to  the 
generation  and  propagation  of  the  sigma  points.  The  number  of  sigma  points  (N)  is  a 
function  of  the  size  of  the  state  (n),  N  =  2n  + 1 ,  and  depending  on  the  size  of  the  state,  the 
computational  load  can  represent  a  limiting  factor. 

2.  Non-additive  Gaussian  Noise 

A  different  model  considering  non-additive  Gaussian  noise  in  the  measurement 

equation  can  be  used,  and  for  this  approach  the  UUV  state  space  representation  becomes 

x(*  +  l)  =  F(*)x(*)+G(*)u(*)  +  o(*) 

i —  _  (5.36) 

z(k )  =  h[k,x(k),co(k)] . 

The  differences  in  the  presented  EKF  and  UKF  algorithms  in  comparison  with  this  case 
are  small  and,  to  avoid  repetition,  are  not  be  presented  in  this  dissertation.  However,  they 
can  be  found  in  [45], 

F.  TUNING 

The  performance  of  KF-based  algorithms  is  heavily  dependent  on  the  tuning 
parameters.  These  parameters  are  often  set  manually  at  great  cost  of  engineering  time 
leading  to,  generally,  non-optimum  choices. 

In  KF  and  EKF  the  tuning  parameters  are  the  covariance  matrices  of  the 
measurement  and  plant  noises  (R  and  Q)  only.  The  UKF  requires,  additionally,  the  tuning 
of  three  scalar  parameters  ( a,  (5 ,  and  k  )  also  known  as  hyper-parameters. 

The  hyper-parameters  a  and  k  are  related  with  the  spreading  of  the  sigma 
points,  and  the  hyper-parameter  /?  is  related  to  the  distribution  of  the  filtered  states. 
Some  heuristics  to  choose  these  tuning  parameters  have  been  proposed  [40],  [41],  but 
there  is  no  consensus  regarding  such  methods.  An  increasing  number  of  authors  have 
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proposed  the  use  of  optimization  techniques  (or  optimizers)  to  select  the  hyper¬ 
parameters  [46],  [47],  [48],  or  the  covariance  matrices  of  the  noises. 

In  this  approach,  the  KF-based  algorithm  is  used  to  construct  the  objective 
function  (or  cost  function)  /(x)  where,  for  the  most  comprehensive  case  (UKF), 

x  =  [«,/?, at, R,Q]  is  subjected  to  the  feasible  set  Q  such  that  xefi.  The  optimization 

technique  will  find  the  best  vector  x  over  all  possibilities  in  Q  that  minimizes  or 
maximizes  the  objective  function. 

The  objective  function  is  not  unique,  and  in  some  circumstances,  the  choice  of  the 
most  suitable  objective  function  can  be  very  challenging.  The  existence  of  discontinuities 
and  multiple  local  extrema  on  the  objective  function  may  be  a  problem  for  the  optimizer. 

Using  the  tracking  data  of  the  KF-based  algorithm,  it  is  possible  to  define  the 
following  objective  functions: 

a)  Residual  measurement  error 

Represents  the  quadratic  measurement  error  weighted  by  the  associated 
covariance  matrix.  It  is  used  when  the  set  of  measurements  is  reliable.  The  optimization 
technique  should  find  a  vector  that  minimizes  the  objective  function 

f(x)  =  'E(z/c-h(^))Sk(^-Hxkjk))  ,  (5.37) 

k=l 

where  zk  is  the  measurement  at  discrete  time  k,  Sk  is  the  covariance  matrix  associated 
with  this  measurement,  h(xk^k )  is  the  nonlinear  measurement  equation  evaluated  for  the 
updated  state  xk^k ,  and  N  is  the  number  of  measurements  available. 

b)  Residual  prediction  error 

Represents  the  quadratic  prediction  error  weighted  by  the  associated  covariance 
matrix.  It  is  used  when  “ground  truth”  data  is  available.  The  optimization  technique 
should  find  a  vector  that  minimizes  the  objective  function 


where  xk  is  the  true  state  at  discrete  time  k,  is  the  updated  predicted  state  at  discrete 

time  k,  P,  is  the  covariance  matrix  associated  with  this  state,  and  N  is  the  number  of  states 
available. 

c)  Log-likelihood  of  the  measurements: 

Used  in  Gaussian  processes  when  the  set  of  measurements  is  reliable.  The 
optimization  technique  should  find  a  vector  that  maximizes  the  objective  function  [49] 

N 

f(x)  =  lnp(zhN\x)  =  YJ 

k=  1 

where  zk  is  the  measurement  at  discrete  time  k,  Sk  is  the  covariance  matrix  associated 
with  this  measurement,  and  N  is  the  number  of  measurements  available. 

d)  Log-likelihood  of  the  states 

Used  in  Gaussian  processes  when  “ground  truth”  data  is  available.  The 
optimization  technique  should  find  a  vector  that  maximizes  the  objective  function  [49] 

N 

f(x)  =  lnp(xhN\x)  =  YJ 

k= 1 

where  xk  is  the  true  state  at  discrete  time  k,  P;  is  the  covariance  matrix  associated  with 
the  predicted  state  x^ ,  and  N  is  the  number  of  states  available. 

In  all  the  preceding  objective  functions,  the  input  data  of  the  KF-based  algorithm 
is  used  as  training  data.  The  quality  of  the  resulting  tuning  will  depend  on  how  well  the 
training  data  represents  the  situations  that  the  filter  could  encounter  in  real  applications. 
Ideally,  an  ample  set  of  training  data,  representing  different  conditions  (ambient  noise, 
sea  state,  relative  positioning  between  assets,  etc.)  should  be  available  to  select  the  tuning 
parameters.  In  environments  where  the  conditions  during  the  mission  are  expected  to 
change  significantly,  the  tuning  would  probably  have  to  be  redone.  In  this  case,  different 
sets  of  tuning  parameters  could  be  stored  in  the  UUV  and  used  accordingly. 

The  selection  of  the  optimal  tuning  parameters,  for  a  given  training  data,  involves 
several  runs  of  the  KF-based  algorithm  (UKF  or  EKF)  and  it  can  take  time  to  converge. 

As  seen  in  Figure  30,  the  process  starts  by  initializing  the  tuning  parameters.  The  UKF 
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It,  1  ,  ,  A’ 

--x[P;1xt--ln|P,|--ln2^ 


(5.40) 


-\zTkS;lzk-hn\Sk  \~\n2n 


(5.39) 


then  runs  and  the  cost  function  is  calculated.  To  be  able  to  find  the  cost  function’s 
extrema  (minimum  or  maximum),  the  optimizer  will  run  UKF  and  calculate  the  cost 
function  for  several  different  tuning  parameters.  Different  techniques  can  be  used  to  find 
the  extrema  and  most  popular  solvers  make  use  of  variational  calculus  algorithms,  as  in 
MATLAB  fmicon  (constrained  nonlinear  minimization),  or  artificial  intelligence 
techniques,  as  in  the  genetic  algorithms  or  simulated  annealing  algorithms. 


Figure  30.  Tuning  Flow  Chart  for  UKF 
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G.  SMOOTHING 


Smoothing  is  a  non-real-time  processing  that  uses  all  measurements  between  ti 
and  ?2  to  estimate  the  state  of  a  system  at  certain  time  t,  where  ty<t<t2  [9],  It  is 

particulary  useful  to  provide  realistic  trajectories  during  the  transient  portions  of  the 
tracking  and  in  situations  where  few  measurements  are  available. 

Smoothing  requires  a  backward  iteration  after  the  forward  filtering  has  been 
performed  resulting  in  x(k|k),  x(k  +  l|k),  P(k|k),  and  P(k  +  l|k)that  should  be  stored 

[38].  The  smoothing  algorithm  presented  in  this  section  is  known  as  the  KF  fixed  interval 
smoother  or  Rauch-Tung-Striebel  smoother. 

Considering  KF  and  EKF,  for  all  k  from  k=N-l  until  k=0  (TV  is  number  of 
samples)  the  following  operations  must  be  performed  [9] : 


a)  Smoother  gain  calculation 

D(k)  =  P(k\k)F(kf  P(k+l\k)~1  . 

b)  Smoothed  state  and  covariance 

x(k|./v)  =  x(k|k)  +  D(k)[^x(k  +  l|./v)-x(k  +  l|k)J 
P(k|A^)  =  P(k|k)  +  D(k)[p(k  +  l|A^)-P(k  +  l|k)]D(k)r  . 


(5.41) 


(5.42) 


For  UKF  additionally,  during  the  forward  filtering,  the  cross-covariance  Ct+1 


must  be  stored  resulting  in  the  following  operations  from  k=N-l  until  k=0  [45]: 


a)  Smoother  gain  calculation 

DW  =  C,.,P(*  +  1|*)_'. 

b)  Smoothed  state  and  covariance 

x(k|  A)  =  x(k|k)  +  D(k)|^x(k  +  l|  Af)-x(k  +  l|k)J 
P(k|A^)  =  P(k|k)  +  D(k)[p(k  +  l|A^)-P(k  +  l|k)]D(k)r  . 


(5.43) 


(5.44) 


H.  CONSENSUS  CURRENT  ALGORITHM 


As  shown  in  Chapter  II,  as  part  of  the  UUV  tracking  model,  the  drift  caused  by 

the  sea  current  was  modeled  as  a  random  walk  and  is  part  of  the  state  of  the  system. 
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Predictions  for  the  sea  current  from  different  UUVs  may  be  transmitted  via  the  reference 
point  at  the  surface  to  a  command  ship  or  shore-based  command  center  for  further 
processing.  The  purpose  of  this  section  is  to  develop  a  consensus  algorithm  able  to  take 
advantage  of  this  information. 

In  the  classical  problem  of  consensus,  a  group  of  agents  need  to  agree  (or  a 
solution  needs  to  converge)  upon  a  certain  quantity  of  interest  [50].  It  is  important  to 
consider  that  some  of  the  agents  may  have  access  to  better  information  than  others.  In  this 
situation,  the  consensus  algorithm  needs  to  be  biased  in  favor  of  those  with  more  reliable 
information  [51]. 

There  are  basically  three  types  of  networks  in  consensus  problems: 

1.  In  a  centralized  network,  the  agents  involved  only  have  the  ability  to 
communicate  with  a  centralized  location; 

2.  In  a  decentralized  network,  the  agents  involved  only  have  the  ability  to 
communicate  with  each  other; 

3.  In  a  mixed  network,  the  agents  have  the  ability  to  communicate  with  a 
centralized  location  and  with  each  other. 

Let  us  expand  the  scenario  described  in  Chapter  II  (depicted  in  Figure  1)  for  a 
multi-agent  (multiple  UUVs)  centralized  network  where  the  assets  at  the  surface  can 
share  information  collected  from  the  UUVs  with  a  shore-based  command  center.  After  a 
certain  time  underwater,  the  UUVs  are  able  to  update  their  position  without  going  to  the 
surface  by  running  the  KF-based  algorithm.  During  this  process,  when  a  measurement  to 
a  surface  asset  is  taken,  the  UUV  may  transmit  its  prediction  for  the  sea  current  and 
associated  uncertainty. 

This  information  is  then  transmitted,  by  the  surface  asset,  to  the  shore  command 
center  using  satellite  communication.  The  command  center  will  receive  the  information 
about  the  sea  current  coming  from  different  UUVs  at  different  times  and  will  need  to 
solve  a  recursive  problem  to  find  the  consensus  sea  current  (see  Figure  31). 
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Iridium  satellite 
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Figure  3 1 .  Centralized  Multi-agent  Network 

The  tracking  algorithm,  running  in  the  UUVs,  is  able  to  make  predictions  for  sea 
current  in  x  and  y  directions,  as  well  as  predict  the  associated  covariance.  We  can  select 
part  of  the  state  and  covariance  related  with  the  sea  current  to  form  the  following 
quantities 

C  =  [C*  S]’  Pe  = 

where 

c  is  the  sea  current  vector  and  Pc  is  the  associated  covariance  matrix; 
er  is  the  covariance  in  x,  er  is  the  covariance  in  y,  and  crtv  =  <jyx  is  the  cross¬ 
covariance. 

Note  that  the  covariance  matrix  is  symmetric.  An  important  property  of 
symmetric  matrices  is  that  the  eigenvalues  ( A )  are  real  and  the  eigenvectors  ( T  )  are 
orthogonal  to  each  other  [52],  Eigenvalues  and  eigenvectors  of  a  covariance  matrix  may 
be  geometrically  interpreted  as 

Eigenvectors  are  orthogonal  unit  vectors  that  represent  the  axis  of  a 
confidence  ellipse  for  which  the  center  is  the  mean  (c); 

Eigenvalues  are  real  numbers  representing  the  length  of  the  ellipse  axes 
(see  Equation  (5.46)  and  Figure  32). 


cr. 


2 


(5.45) 
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PCY  =  W 


e  0 

COS# 

-sin# 

=>  A  = 

i 

o 

i 

;  Y  = 

sin# 

cos# 

(5.46) 


sin0) 


( is  a  scaling  factor  to  attribute  a  certain  percentage  to  the  confidence  ellipse,  for  example, 
for  a  95%  confidence  ellipse  (considering  normally  distributed  data)  (=  V5.991. 

Figure  32.  Ellipse  Representing  a  Confidence  Interval  around  the  Mean  c 

According  to  Equation  (5.46) 


We  can  say  that  the  eigenvector  associated  with  the  smaller  eigenvalue  points  to 
the  direction  with  smaller  variation  around  the  mean  (or  points  to  a  more  reliable  region). 
From  Equation  (5.46)  and  Figure  32,  the  eigenvector  associated  with  the  smaller 
eigenvalue  is 


H  = 


-sin# 

cos# 


(5.47) 


The  eigenvector  H  may  be  used  to  reduce  the  amount  of  data  transmitted  by  the  UUV, 
while  still  carrying  information  about  the  uncertainty  associated  with  the  prediction  for 
the  sea  current.  This  may  be  done  by  writing 


y=cH  =  [c, 


1  -sin# 
cos# 


(5.48) 


The  UUV  may  now  transmit  only  two  real  numbers  (the  angle  #  and  the  quantity 
y)  instead  of  five  real  numbers  as  in  Equation  (5.45).  The  command  center  will  receive  y 
and  H  from  different  UUVs  at  different  times  and  have  to  solve  recursively  for  c  ( k ) , 


using  the  relationship 
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y(k)  =  c(k)H(k).  (5.49) 

Equation  (5.49)  may  be  thought  of  as  being  a  measurement  and  the  following  state  space 
representation  may  be  used 


(5.50) 


c(k  +  l)  =  F(k)c(k)  +  v(k) 
y(k)  =  H(k)c(k)  +  oi(k), 

where  F  ( k )  is  the  identity  matrix  and  the  noises  o ( k )  and  co(k)  are  assumed  to  be  zero 

mean,  white  (i.e.,  uncorrelated)  Gaussian  processes. 


The  sea  current  c(k)  in  Equation  (5.50)  may  be  solved  using  the  Kalman 

equations  as  shown  in  Section  A.  When  the  predictions  for  the  current  converge  to  a 
steady-state  value,  we  have  the  consensus  sea  current  (CSC).  Uses  for  the  CSC  may  be  as 
follows: 

the  command  center  can  use  CSC  in  its  upper  level  tracking  algorithm; 
CSC  can  be  broadcasted  to  all  UUVs  to  be  used  in  their  tracking 
algorithms; 

CSC  can  help  to  decide  if  two  UUVs  have  reliable  position  to  take 
distance  measurements  from  each  other  (CSC  amplitude  and  time  of  the 
last  position  update  could  be  the  criteria); 

A  complete  view  of  the  consensus  algorithm  is  shown  in  Figure  33. 
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Figure  33.  Consensus  Current  Algorithm 
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VI.  SIMULATIONS 


In  this  chapter,  the  algorithms  developed  in  Chapter  V  are  tested  considering  three 
distinct  cases.  The  first  is  the  ideal  case  where  all  the  measurements  are  successful  and  no 
noise  is  present.  The  second  case  considers  failure  in  the  measurements  during  a  certain 
portion  of  the  UUVs  mission.  In  the  third,  noise  is  added  to  the  previous  case. 

A  scenario  in  which  two  UUVs  navigating  in  an  area  where  three  surface  assets 
are  present  is  simulated.  The  surface  assets,  represented  by  Wave  Gliders  (WG)  in 
Figure  34,  have  access  to  GPS  and  can  share  information  with  the  shore-based  command 
center  via  satellite  as  in  a  centralized  multi-agent  network,  permitting  evaluation  of  the 
consensus  current  algorithm. 


UUV-2 


Figure  34.  Simulation  Scenario 

The  true  trajectory  of  the  UUVs  is  depicted  in  Figure  35.  The  UUVs  are  set  to 
navigate  at  1  knot  in  an  area  with  a  constant  current  of  0.5  knots  at  135°.  The  starting 
point  for  UUV-1  is  (-400m,  1400m)  and  for  UUV-2  is  (1000m,  1400m),  in  a  Cartesian 
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coordinate  system.  For  simplicity,  the  WGs  are  considered  stationed  at  fixed  positions,  as 
in  a  perfect  hovering  mission. 


Figure  35.  UUVs  True  Trajectory 


A.  IDEAL  CASE 

In  this  ideal  case,  each  UUV  has  a  successful  measurement  every  three  minutes, 
as  depicted  in  Figure  36.  As  described  in  Chapter  II,  a  measurement  from  the  UUVs  to 
two  or  more  reference  points  (WGs)  at  the  same  time  is  not  possible,  and  no  bearing 
information  associated  with  the  distance  measurements  is  available.  Therefore,  in  the  xy- 
plane,  each  measurement  represents  a  circle  of  possible  UUV  positions  (see  Figure  37). 
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Figure  36.  Time  Evolution  of  the  Distance  Measurements,  Ideal  Case 


85 


-1500  -1000  -500  0  500  1000  1500  2000  2500  3000 

Xaxis  [m] 


Xaxis  [m] 


Figure  37.  Distance  Measurements  in  xy-Plane,  Ideal  Case 

Figure  38  shows  the  tracking  results  (predicted  and  smoothed  trajectories)  for 
both  algorithms  (EKF  and  UKF)  and  in  Figure  39  the  tracking  errors  are  presented.  The 
prediction  errors  in  Figure  39  show  that  both  algorithms  converge  to  the  true  trajectory, 
but  UKF  converges  faster.  A  comparison  in  Figure  39  between  predictions  and  smoothed 
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trajectories  shows  how  useful  the  smoother  is  to  reduce  the  errors  during  the  transient 
portion  of  the  tracking. 

Figure  40  shows  the  consensus  current  results.  As  UKF  running  in  the  UUVs  was 
able  to  converge  to  the  true  current  faster,  the  consensus  current  based  on  the  UKF 
tracking  has  a  smoother  profile  and  converges  faster  than  EKF. 


Figure  38.  Tracking  Results,  Ideal  Case 
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UUV-1  Trajectory  Errors 


UUV-2  Trajectory  Errors 


UUV-1  Trajectory  Errors 


UUV-2  Trajectory  Errors 


Figure  39.  Error  Analysis,  Ideal  Case 
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Figure  40.  Consensus  Current,  Ideal  Case 


B.  MEASUREMENT  FAILURE 

In  this  case,  the  UUVs  are  experiencing  a  failure  in  the  distance  measurement 
during  part  of  their  mission,  as  depicted  in  Figures  41  and  42.  The  tracking  results  show 
that  both  algorithms  converge  to  the  true  trajectory,  although  the  EKF  takes  more  time 
and  is  highly  affected  during  the  turn  (see  Figures  43  and  44).  A  comparison  in  Figure  44 
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between  predictions  and  smoothed  trajectories  shows  that  the  smoother  was  able  to 
reduce  the  tracking  errors  during  the  transient  portion  of  the  tracking  and  during  the 
period  of  failure  in  the  measurements. 

Figure  45  shows  the  consensus  current  results.  As  in  the  previous  case,  because 
UKF  running  on  the  UUVs  was  able  to  converge  to  the  true  current  faster,  the  consensus 
current  based  on  UKF  tracking  has  a  smoother  profile  and  converges  faster  than  EKF. 


Figure  41.  Time  Evolution  of  the  Distance  Measurements, 

Measurement  Failure  Case 
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Figure  42.  Distance  Measurements  in  xy-Plane,  Measurement  Failure  Case 
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Figure  43.  Tracking  Results,  Measurement  Failure  Case 


92 


UUV-1  Trajectory  Errors  UUV-1  Trajectory  Errors 


UUV-2  Trajectory  Errors  UUV-2  Trajectory  Errors 


Figure  44.  Error  Analysis,  Measurement  Failure  Case 
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EKF  Currents 
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Figure  45.  Consensus  Current,  Measurement  Failure  Case 


C.  SIMULATIONS  WITH  MEASUREMENT  NOISE 

For  this  case,  in  addition  to  measurement  failures  presented  in  Figures  41  and  42, 
noise  is  added  in  the  distance  measurements.  To  account  for  the  random  characteristics  of 
the  noise,  Figures  46,  47,  and  48  are  the  result  of  the  average  of  100  simulations. 
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As  described  in  Chapter  IV,  distances  underwater  are  estimated  using  the  one-way 
travel  time,  and  the  one-way  travel  time  is  estimated  by  the  time  of  the  highest  peak  in 
the  modem’s  matched  filter  output.  In  a  multipath  environment,  the  arrival  with  higher 
energy  (highest  peak)  is  not  always  the  direct  path  (first  group  of  arrivals)  but  is 
representative  of  a  bottom/surface  bounce  arriving  later.  As  a  consequence,  the  distance 
is  overestimated,  which  means  that  the  errors  associated  with  the  distance  measurements 
are  always  positive  (biased  noise).  This  particular  characteristic  is  considered  in  the 
simulations. 

The  tracking  results  show  that  both  algorithms  converge  to  the  true  trajectory  but 
the  errors  are  higher  than  in  the  previous  cases.  In  Figure  47,  it  can  be  seen  that  UKF  has 
a  more  stable  behavior  and  smaller  errors  when  compared  with  EKF,  probably  because 
UKF  was  able  to  handle  the  biased  noise  better.  This  result  is  consistent  with  the  claim 
that  UKF,  by  the  use  of  sigma  points,  can  capture  the  characteristics  of  the  true 
distribution  better  than  EKF. 

A  comparison  in  Figure  47  between  predictions  and  smoothed  trajectories  shows 
that  the  smoother  was  able  to  reduce  the  tracking  errors.  In  Figure  48,  the  consensus 
current  results  show  that  the  consensus  based  on  the  UKF  tracking  has  a  smoother  profile 
and  converges  faster  than  EKF. 
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Figure  46.  Tracking  Results,  Measurement  Failure  Case  Plus  Noise 
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UUV-1  Trajectory  Errors  UUV-1  Trajectory  Errors 


UUV-2  Trajectory  Errors  UUV-2  Trajectory  Errors 


Figure  47.  Error  Analysis,  Measurement  Failure  Case  Plus  Noise 
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EKF  Currents 
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Figure  48.  Consensus  Current,  Measurement  Failure  Case  Plus  Noise 

D.  CURRENT  ESTIMATION  BY  THE  CONSENSUS  ALGORITHM 

In  the  previous  sections,  the  consensus  current  algorithm  was  demonstrated  to 
produce  a  smoother  profile  and  quicker  convergence  than  EKF  when  the  tracking  results 
from  the  UKF-based  algorithm  were  used  as  input.  In  this  section,  the  tracking  results 
from  UKF  are  used  in  a  situation  where  the  UUVs  are  experiencing  slightly  different 
currents.  UUV-1,  navigating  at  1  knot,  is  experiencing  a  current  of  0.45  knots  at  135°  and 
UUV-2,  navigating  at  1  knot,  is  experiencing  a  current  of  0.6  knots  at  132°.  Their 
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trajectories  are  depicted  in  Figure  49.  The  scenario  defined  was  consistent  with  the  case 
described  in  Section  C  (measurement  failure  plus  noise). 
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Figure  49.  True  Trajectory  of  UUVs  (Experiencing  Different  Currents) 


UKF  Currents 


Figure  50.  Consensus  Current,  UUVs  Experiencing  Different  Currents 


In  Figure  50  it  can  be  seen  that  the  consensus  current  stabilizes  to  a  value  between 
the  two  true  currents.  Probably  due  to  a  lower  covariance  associated  with  the  current 
predictions  from  UUV-2  (i.e.,  predictions  for  UUV-2  converge  faster),  the  consensus 
current  is  closer  to  this  value. 
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VII.  SEA  TESTS 


In  this  chapter,  the  results  of  the  sea  tests  that  took  place  in  Monterey  Bay  in 
August  2015  are  presented.  In  a  one-week  sea  test,  two  Liquid  Robotics  Wave  Gliders 
(unmanned  surface  vehicles,  or  US  Vs  named  Mako  and  Tiburon)  and  a  command  ship 
(NOOA  R/V  Fulmar),  all  equipped  with  Teledyne-Benthos  acoustic  modems  (ATM-900 
series)  and  GPS  connectivity  are  forming  a  network  of  reference  points  for  an  Exocetus 
Coastal  Glider  (UUV  named  LG  16),  which  is  also  equipped  with  the  same  type  of 
acoustic  modem  (see  Figure  51). 

The  objective  of  the  sea  test  was  to  evaluate  the  ability  of  the  tools  developed  in 
this  work  (i.e.,  the  tracking  algorithm  and  matched  peak  impulse  response  processor)  to 
track  the  UUV.  During  the  sea  tests,  several  one-hour  missions  were  successfully 
conducted.  For  each  mission,  the  UUV  navigation  system  recorded  its  position 
predictions  and,  due  to  miscalibration  and  misalignment  of  its  sensors  and  the  sea  current, 
errors  of  the  order  of  500  meters  and  above  were  observed  between  the  final  UUV 
predicted  position  and  its  GPS  location,  measured  soon  after  it  surfaced. 

During  the  course  of  the  missions,  several  travel  time  measurements  between  the 
UUV  and  the  reference  points  were  successfully  recorded.  The  travel  times  were 
converted  to  horizontal  distances  and  supplied  to  the  tracking  algorithm  for  processing. 


Figure  5 1 .  Sea  Test  Configuration 
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A. 


KEY  ASSETS 


This  section  provides  the  basic  characteristics  and  capabilities  of  the  UUV  and  the 
Wave  Gliders  used  in  the  sea  trials. 

1.  Exocetus  Coastal  Glider 

This  vehicle  makes  use  of  a  buoyancy  engine  to  achieve  forward  motion  in  a  saw 
tooth  pattern  by  descending  and  ascending  from  changes  in  its  buoyancy.  There  are  no 
external  moving  parts  or  propeller.  To  maneuver,  the  glider  has  to  pitch  and  roll  its 
internal  battery.  It  can  carry  a  payload  of  5  kg,  and  can  reach  speeds  up  to  2  knots.  Its 
navigation  system  relies  on  a  digital  magnetic  compass  board  and  speed  estimates  based 
on  the  buoyancy  engine  dynamics.  Communications  capabilities  include  Iridium,  UHF 
(Freewave  radio  modem),  and  WiFi. 


Figure  52.  Exocetus  Coastal  Glider.  Source:  Exocetus  Autonomous  Systems, 

www.exocetussystems.com 


Here  at  the  Naval  Postgraduate  School  the  vehicle  was  rigged  on  its  head  with  a 
syntactic  foam  mount  for  the  acoustic  modem  as  seen  in  Figure  53.  The  purpose  of  the 
synthetic  foam  mount  was  to  offset  the  weight  of  the  modem  transducer  to  avoid 
interference  with  the  buoyancy  engine  operations. 
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Figure  53.  Coastal  Glider  Rigged  with  Syntactic  Foam  Mount  for  Acoustic  Modem 

2.  Liquid  Robotics  Wave  Glider 

This  is  a  hybrid  vehicle  composed  of  a  surface  (or  float)  unit  and  an  underwater 
(or  sub)  unit  attached  to  each  other  via  a  tether  (umbilical  cable).  The  underwater  unit  is 
able  to  convert  part  of  the  upward  and  downward  motion  of  the  water  column  into 
forward  motion  by  using  its  cantilevered  wings.  As  waves  move  the  system  up  and  down, 
the  underwater  unit  acts  as  a  tug  pulling  the  surface  float  along  a  predetermined  course.  It 
can  tow  up  to  500  kg,  and  can  reach  speeds  up  to  3  knots.  It  has  GPS  connectivity  and 
communications  capabilities  that  include  Iridium  and  WiFi.  The  WGs  used  in  the  sea 
tests  also  have  an  acoustic  modem  installed  in  a  towfish  structure,  for  communication 
with  the  UUV,  as  represented  schematically  in  Figure  55. 


Figure  54.  Liquid  Robotics  Wave  Glider  Model  SV2.  Source:  Liquid  Robotics, 

https://www.liquid-robotics.com/ 


103 


T  owfish 


Acoustic  Modem 

Figure  55.  Wave  Glider  Equipped  with  Towfish  Modem.  Adapted  from  Liquid  Robotics, 

https://www.liquid-robotics.com/ 

B.  AREA  OF  OPERATION  AND  ENVIRONMENTAL  CHARACTERISTICS 

The  assets  were  operating  in  a  relatively  flat  area  around  75  m  of  depth  at  the 
edge  of  the  Monterey  Bay  submarine  canyon  about  5  nautical  miles  from  the  shore  (see 
Figure  19).  Sediments  in  the  area  were  identified  as  being  silt-clay/sand-silt-clay  and 
modeled  as  a  fluid-like  half-space,  according  to  [37],  with  the  following  characteristics: 

Compressional  sound  speed:  1560  m/s; 

Density:  1.6  g/cm3; 

Attenuation  for  compressional  waves:  5dB/m. 

The  sound  speed  profile  was  measured  by  the  UUV  during  the  tests,  and  its  average  is 
depicted  in  Figure  56. 
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Figure  56.  Average  Sound  Speed  Profile  Measured  by  the  UUV  on  August  12,  2015 

C.  UUV  MISSIONS 

After  a  week  of  preliminary  tests,  on  August  12,  2015,  four  UUV  missions  were 
run.  For  three  of  them,  successful  travel  time  measurements  from  the  UUV  to  the 
reference  points  were  recorded.  Those  travel  times  were  converted  into  distances  using 
the  algorithm  described  in  Chapter  IV,  Section  B.2.c.  The  travel  time  data,  the  data  from 
the  UUV  navigation  system,  and  the  coordinates  of  the  reference  points  were  supplied  to 
the  tracking  algorithm  for  processing.  All  data  collected  was  post  processed. 

The  tuning  parameters  were  selected  automatically  using  the  algorithm  described 
in  Chapter  V,  Section  E,  and  the  objective  function  used  was  the  residual  prediction  error. 
Data  after  the  smoothing  was  used  to  calculate  the  objective  function.  Unfortunately,  the 
EKF  method  proved  to  be  unstable  for  automatic  tuning  using  real  data.  The  results 
presented  in  this  section  include  only  UKF. 

To  verify  the  effectiveness  of  the  algorithm  described  in  Chapter  IV,  Section 
B.2.c,  a  comparison  of  the  tracking  results  using  distance  estimations  according  to  that 
section  and  distances  roughly  estimated  are  performed.  For  the  purpose  of  this 
dissertation,  the  term  distances  roughly  estimated  is  defined  as  being  the  result  of  the 
multiplication  of  the  one-way  travel  time  by  the  characteristic  medium  sound  speed, 
corrected  by  the  asset’s  depth  (representing  a  horizontal  distance). 
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Based  on  the  sound  speed  profile  depicted  in  Figure  56,  the  rough  estimates  of  the 
distances  from  the  UUV  to  the  reference  points  are  calculated  using  a  constant  sound 
speed  of  1465  m/s  (corresponding  to  the  sound  speed  associated  with  typical  depths  of 
the  reference  points  and  the  UUV  when  the  measurements  were  recorded). 

1.  Data  Collection  and  Sources  of  Errors 

Despite  our  efforts  to  minimize  and  account  for  errors  in  measurements,  we 
identified  the  following  main  sources  of  errors  that  introduce  uncertainties  and  may  affect 
the  tracking  results: 

-  Wave  Glider  GPS  positions:  The  WGs  record  their  GPS  positions  in  5 -minute 
intervals.  In  order  to  establish  their  position  at  the  instant  of  a  travel  time 
measurement  it  becomes  necessary  to  interpolate  between  GPS  recordings. 

-  Offset  between  the  GPS  antenna  and  the  acoustic  modem:  Due  to  the  use  of  the  tow 
fish,  the  GPS  antenna  on  a  WG  is  displaced  from  the  acoustic  modem  transducer 
(see  Figure  55).  The  estimated  displacement  is  around  4  to  6  meters,  and  the  sub 
unit  heading  can  be  used  to  apply  a  correction  in  the  proper  direction. 

-  Sound  speed  profile:  The  sound  speed  profile  represented  in  Figure  56  is  the 
average  of  the  UUV’s  measurements  during  several  dives  throughout  its  mission 
(and  extrapolated  to  the  area  depth).  The  actual  measurements  are  depicted  in 
Figure  57.  The  variability  in  the  sound  speed  profile  may  cause  errors  in  distance 
estimations  based  on  the  matched  peak  IR  algorithm  by  affecting  the  amplitude  of 
the  arrivals. 

-  Acoustic  modem  beam  patterns:  As  mentioned  in  Chapter  IV,  Section  B.2.d,  the 
vertical  and  horizontal  beam  patterns  are  considered  omnidirectional.  To  be  able  to 
take  into  account  the  standard  transducer  beam  pattern,  as  well  as  the  effects  of 
mounting  onto  the  UUV  and  tow  fish  bodies,  the  beam  patterns  should  be  measured 
when  the  modems  are  mounted  in  their  respective  structures.  The  omnidirectional 
assumption  may  cause  errors  in  distance  estimations  based  on  the  matched  peak  IR 
algorithm  by  affecting  the  amplitude  of  the  arrivals. 
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Figure  57.  Sound  Speed  Profile  as  Measured  by  the  UUV  on  August  12,  2015 

2.  Mission-1 

Mission-1  started  at  09:34:24  AM  and  ended  at  10:36:48  AM.  During  this  period 
of  time  the  UUV  navigation  system  recorded  its  estimated  postion.  Due  to  errors  in  dead 
reckoning  and  the  sea  current,  an  error  of  1095  meters  was  observed  between  the  final 
UUV  predicted  position  and  its  GPS  location  measured  soon  after  it  surfaced  (see  Figure 
58). 

During  the  course  of  this  mission,  several  travel  time  measurements  between  the 
UUV  and  the  reference  points  were  successfully  recorded.  The  travel  times  were 
converted  to  horizontal  distances  using  the  algorithm  described  in  Chapter  IV,  and  the 
evolution  of  the  distance  measurements  over  time  is  depicted  in  Figure  59. 

Note  that  the  measurements  were  recorded  from  the  UUV  to  only  two  reference 
points  and  in  an  inconsistent  manner.  Ideally,  the  measurements  should  be  taken  at 
constant  time  intervals  to  different  reference  points  in  a  orderly  manner,  as  shown  in 
Chapter  VI,  Figure  36.  Another  important  point  is  that  no  measurement  was  recorded 
during  the  first  3  minutes  and  30  seconds,  and  during  the  last  5  minutes  and  43  seconds  of 
the  mission,  as  well. 
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08/12/2015  Tests,  Mission-1  -  Origin  Lat:36.79  Lon:-121.89 


Figure  58.  August  12,  2015,  Mission- 1 


Time 


Figure  59.  Mission- 1,  Time  Evolution  of  the  Distance  Measurements 


The  tracking  results  after  automatic  tuning,  using  the  distances  from  the  UUV  to 
the  reference  points  estimated  according  Chapter  IV,  Section  B.2.c,  are  presented  in 
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Figure  60.  Despite  the  inconsistent  measurements,  the  automatic  tuning  algorithm  was 
able  to  find  the  optimum  parameters  that  provided  consistent  tracking. 


08/12/2015  Tests,  Mission-1  -  Origin  Lat:36.79  Lon:-121.89 


Figure  60.  Mission- 1,  Results  after  Automatic  Tuning 
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As  seen  in  Chapter  V,  Section  G,  for  each  prediction  there  is  an  associated 
covariance  matrix,  and  the  eigenvalues  and  eigenvectors  of  the  covariance  matrices  may 
be  used  to  represent  confidence  ellipses  around  the  predictions.  Despite  some  limitations 
[38],  this  procedure  can  provide  good  information  about  the  uncertainty  associated  with 
the  predicted  states. 

In  Figure  61,  the  smoothed  covariance  matrix  associated  with  the  last  prediction 
was  used  to  construct  a  confidence  ellipse  around  it.  The  lengths  of  the  main  axes  were 
calculated  by  the  square  root  of  the  associated  covariance  matrix  eigenvalues  and  the 
main  axes’  directions  were  taken  from  the  eigenvectors.  The  confidence  ellipse  major 
axis  is  7  meters,  and  the  minor  is  2  meters.  Despite  the  numbers  in  Table  2  indicating  a 
“perfect  tracking”  (i.e.,  the  last  prediction  matches  the  GPS  fix),  the  most  important  result 
shown  by  Figure  61  is  the  realistic  size  of  the  confidence  ellipse  and  the  fact  that  the  GPS 
fix  is  located  inside  this  ellipse. 
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Figure  61.  Mission- 1,  Confidence  Ellipse  at  the  End  of  the  Mission 


The  optimum  tuning  parameters  for  this  mission  are  shown  in  Table  1.  The 

parameter  *J~R  represents  the  uncertainty  associated  with  the  measurements,  combining 
the  uncertainties  in  the  horizontal  distance  and  in  the  position  of  the  reference  points  (in  x 
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and  y).  A  large  *J~R  may  represent  a  non-reliable  or,  in  this  case,  inconsistent 
measurement.  As  already  pointed  out,  in  Figure  59  can  be  seen  that  we  have  scattered 
measurements  to  only  two  of  the  three  reference  points. 


Table  1.  Mission- 1,  Optimum  Tuning  Parameters  for  Input  Data  Estimated 

According  to  Chapter  IV,  Section  B.2.d 


Tuning  Parameters 

a 

1.08X10'1 

p 

3.82 

K 

2.03  xlO"9 

q  [m2/s3]a 

2.09  xlO"3 

VR  [m2] 

4.65xl05 

a,  P  ,  and  K  have  no  dimensions; 
a  Units  according  to  [38], 


For  this  mission,  tracking  by  dead  reckoning  (UUV  navigation  system)  yielded  an 
error  of  1095  meters  (see  Figure  58).  The  tracking  algorithm  was  tested  and  the  results 
show  a  more  reliable  predicted  position  for  the  entire  UUV  trajectory  (see  Figure  60). 

When  a  rough  estimate  for  the  distances  is  used,  errors  around  30  meters  were 
observed  at  the  end  of  the  mission.  While  this  is  a  significant  improvement  over  the 
UUV’s  estimated  position,  it  is  based  on  simple  straight  line  propagation  that  does  not 
account  for  the  refractive  and  multipath  effects  of  the  environment  (Figure  62). 

When  the  more  complete  algorithm  described  in  Chapter  IV  is  used,  the  tracking 
algorithm  practically  eliminates  the  error  at  the  end  of  the  mission  (as  represented  in 
Figure  61,  there  are  uncertainties  associated  with  this  prediction).  This  indicates  that 
correlating  the  measured  impulse  response  of  the  channel  with  a  ray  trace  prediction  may 
provide  an  improvement  in  the  tracking  accuracy. 
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Table  2.  Mission- 1,  Tracking  Comparison 


Mission- 1 

Distances  according 
to  Chapter  IV 

Distances 
roughly  estimated 

Objective 

Function  [m4] 

0.47 

3.01xl08 

Error  at  beginning 
of  the  Mission  [m] 

0.04 

5.20 

Error  at  end 
of  the  Mission  [m] 

0 

29.7 

Total  Error  [m] 

0.04 

34.9 

The  row  Error  at  beginning  of  the  Mission  on  Tables  2,  4,  and  6  refers  to  the 
difference  between  the  first  predicted  position  and  the  UUV  GPS  fix  at  the  beginning  of 
the  mission.  The  row  Error  at  end  of  the  Mission  refers  to  the  difference  between  the  last 
predicted  position  and  the  UUV  GPS  fix  at  the  end  of  the  mission. 
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Figure  62.  Mission- 1,  Tracking  Results  when  Processing  the  Distances  Estimated  According 
to  Chapter  IV  and  when  Processing  the  Distances  Roughly  Estimated 

The  sea  current  predicted  by  the  tracking  algorithm  for  this  mission  is  0.60  knots 
(i.e.,  the  average  of  the  smoothed  predictions)  and  has  a  general  direction  as  represented 
in  Figure  63a.  It  is  worth  noting  the  evolution  of  the  tide  during  the  mission  period  (see 
Figure  63b). 
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Figure  63.  Mission- 1,  Tracking,  Sea  Current,  and  Tide 


3.  Mission-2 

Mission-2  started  at  10:49:03  AM  and  ended  at  11:51:38  AM.  During  this  period 
of  time  the  UUV  navigation  system  recorded  its  estimated  position.  Due  to  errors  in  dead 
reckoning  and  the  sea  current,  an  error  of  747  meters  was  observed  between  the  final 
UUV  predicted  position  and  its  GPS  location,  measured  soon  after  it  surfaced  (see  Figure 
64). 

During  the  course  of  this  mission,  travel  time  measurements  between  the  UUV 
and  the  reference  points  were  successfully  recorded.  The  travel  times  were  converted  to 
horizontal  distances  using  the  algorithm  described  in  Chapter  IV,  and  the  evolution  of  the 
distance  measurements  over  time  is  depicted  in  Figure  65. 

Note,  in  Figure  65,  that  just  a  few  measurements  were  recorded,  but  in  a  more 
consistent  way  when  compared  with  mission- 1.  Another  interesting  point  is  that 
measurements  to  all  three  reference  points  were  successfully  recorded  and  the  last  one 
was  taken  just  2  minutes  and  21  seconds  before  the  UUV  GPS  fix  at  the  end  of  the 
mission.  Note,  too,  that  no  measurement  was  recorded  during  the  first  33  minutes  and  50 
seconds  of  this  mission. 
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Figure  64.  August  12,  2015,  Mission-2 
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Figure  65.  Mission-2,  Time  Evolution  of  the  Distance  Measurements 


The  tracking  results  after  automatic  tuning,  using  the  distances  from  UUV  to  the 
reference  points  estimated,  as  described  in  Chapter  IV,  Section  B.2.c,  are  presented  in 
Figure  66.  Despite  the  few  measurements,  the  automatic  tuning  was  able  to  find  the 
optimum  parameters  that  provided  a  consistent  tracking. 
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Due  to  the  long  period  without  measurements  (at  the  beginning  of  the  mission) 
and  due  to  the  poor  performance  of  the  UUV  navigation  system  predicting  its  velocity, 
the  non-smoothed  predictions  (represented  by  the  black  curve  in  Figure  66b)  appear  to  be 
significantly  off  track  before  the  measurements  (represented  by  the  “jumps”  on  the  black 
curve).  But  measurements  to  different  reference  points  were  a  decisive  factor  in  getting 
the  track  on  the  correct  path.  In  this  particular  case,  the  smoothing  algorithm  proved  very 
important  in  producing  a  realistic  UUV  trajectory. 


1000  1500  2000  2500 

x  [m] 
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Figure  66.  Mission-2,  Results  after  Automatic  Tuning 
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In  Figure  67,  the  smoothed  covariance  matrix  associated  with  the  last  prediction 
was  used  to  construct  a  confidence  ellipse  around  it.  The  lengths  of  the  main  axes  were 
calculated  by  the  square  root  of  the  associated  covariance  matrix  eigenvalues,  and  the 
main  axes’  directions  were  taken  from  the  eigenvectors.  The  confidence  ellipse  major 
axis  is  8  meters,  and  the  minor  is  2  meters.  Despite  the  numbers  in  Table  4  indicating  a 
“perfect  tracking”  (i.e.,  the  last  prediction  matches  the  GPS  fix),  the  most  important  result 
shown  in  Figure  67  is  the  realistic  size  of  the  confidence  ellipse  and  the  fact  that  the  GPS 
fix  is  located  inside  this  ellipse. 
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Figure  67.  Mission-2,  Confidence  Ellipse  at  the  End  of  the  Mission 


The  optimum  tuning  parameters  for  this  mission  are  shown  in  Table  3.  As  can  be 

seen,  probably  due  to  the  more  consistent  measurement,  the  parameter  \[r  is  in 
reasonable  levels. 
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Table  3.  Mission-2,  Optimum  Tuning  Parameters  for  Input  Data  Estimated 

According  Chapter  IV,  Section  B.2.d 


Tuning  Parameters 

a 

3.05x10'* 

fi 

8. 98x1  O'2 

K 

9.40x1 0"3 

q  [m2/s3]a 

l.llxlO'3 

Vr  [m2] 

10 

a,  (3  ,  and  K  have  no  dimensions; 
a  Units  according  to  [38], 


For  this  mission,  tracking  by  dead  reckoning  (UUV  navigation  system)  yielded  an 
error  of  747  meters  (see  Figure  64).  The  tracking  algorithm  was  tested  and  the  results 
show  a  more  reliable  predicted  position  for  the  entire  UUV  trajectory  (see  Figure  66). 
When  a  rough  estimate  for  the  distances  is  used,  errors  around  39  meters  were  observed 
at  the  end  of  the  mission.  While  this  is  a  significant  improvement  over  the  UUV’s 
estimated  position,  it  is  based  on  simple  straight  line  propagation  that  does  not  account 
for  the  refractive  and  multipath  effects  of  the  environment  (Figure  68). 

When  the  more  complete  algorithm  described  in  Chapter  IV  is  used,  the  tracking 
algorithm  produces  no  apparent  error  at  the  end  of  the  mission  (as  represented  in  Figure 
67,  there  are  uncertainties  associated  with  this  prediction).  This  indicates  that  correlating 
the  measured  impulse  response  of  the  channel  with  a  ray  trace  prediction  may  provide  an 
improvement  in  the  tracking  accuracy. 


Table  4.  Mission-2,  Tracking  Comparison 


Mission-2 

Distances  according  to 
Chapter  IV 

Distances 
roughly  estimated 

Objective 

Function  [m4] 

0.04 

402 

Error  at  beginning 
of  the  Mission  [m] 

0.009 

1.02 

Error  at  end 
of  the  Mission  [m] 

0.004 

39.1 

Total  Error  [m] 

0.013 

40.1 
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Figure  68.  Mission-2,  Tracking  Results  when  Processing  the  Distances  Estimated 
According  to  Chapter  IV  and  when  Processing  the  Distances 
Roughly  Estimated 


The  sea  current  predicted  by  the  tracking  algorithm  for  this  mission  is  0.45  knots 
(i.e.,  the  average  of  the  smoothed  predictions)  and  has  a  general  direction  as  represented 
in  Figure  69a.  It  is  worth  noting  the  evolution  of  the  tide  during  the  mission  period  (see 
Figure  69b). 
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Figure  69.  Mission-2  Tracking,  Sea  Current,  and  Tide 
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4. 


Mission-3 


Mission-3  started  at  12:05:01  PM  and  ended  at  13:07:38  PM.  During  this  period 
of  time  the  UUV  navigation  system  recorded  its  estimated  postion.  Due  to  errors  in  dead 
reckoning  and  the  sea  current,  an  error  of  476  meters  was  observed  between  the  final 
UUV  predicted  position  and  its  GPS  location,  measured  soon  after  it  surfaced  (see  Figure 
70). 

During  the  course  of  this  mission,  travel  time  measurements  between  the  UUV 
and  the  reference  points  were  successfully  recorded.  The  travel  times  were  converted  to 
horizontal  distances  using  the  algorithm  described  in  Chapter  IV,  and  the  evolution  of  the 
distance  measurements  over  time  is  depicted  in  Figure  71. 

Note  the  inconsistent  measurement  at  the  beginning  and  the  lack  of  measurements 
just  after  the  middle  of  the  mission.  In  the  last  portion,  however,  measurements  to  all 
three  reference  points  were  achieved  in  a  short  period  of  time.  The  first  measurement  was 
not  taken  until  6  minutes  and  15  seconds  into  the  mission  and  the  last  one  was  taken  6 
minutes  and  37  seconds  before  the  UUV  GPS  fix  at  the  end. 
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Figure  70.  August  12,  2015,  Mission-3 
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Figure  71.  Mission-3,  Time  Evolution  of  the  Distance  Measurements 

The  tracking  results  after  automatic  tuning,  using  the  distances  from  UUV  to  the 
reference  points  estimated  as  decribed  in  Chapter  IV,  Section  B.2.c,  are  presented  in 
Figure  72.  Despite  the  reduced  number  of  measurements  and  the  inconsistency  during 
certain  portions  of  the  mission,  the  automatic  tuning  was  able  to  find  the  optimum 
parameters  that  produced  an  error  of  8.6  meters  at  the  end  of  the  mission  (note  that  there 
are  uncertainties  associated  with  this  prediction,  as  shown  in  Figure  73). 

In  Figure  73,  the  smoothed  covariance  matrix  associated  with  the  last  prediction 
was  used  to  construct  a  confidence  ellipse  around  it.  The  lengths  of  the  main  axes  were 
calculated  by  the  square  root  of  the  associated  covariance  matrix  eigenvalues,  and  the 
main  axes’  directions  were  taken  from  the  eigenvectors.  The  confidence  ellipse  major 
axis  is  27  meters,  the  minor  is  14.5  meters,  and  the  GPS  fix  is  located  inside  this  ellipse. 
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Figure  72.  Mission-3,  Results  after  Automatic  Tuning 
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Figure  73.  Mission-3,  Confidence  Ellipse  at  the  End  of  the  Mission 


The  optimum  tuning  parameters  for  this  mission  are  shown  in  Table  5.  As  can  be 

seen,  probably  due  to  the  more  consistent  measurement,  the  parameter  *J~R  is  in 
reasonable  levels. 


Table  5.  Mission-3,  Optimum  Tuning  Parameters  for  Input  Data  Estimated 
According  to  Chapter  IV,  Section  B.2.d 


Tuning  ] 

>arameters 

a 

1.42x1  O'2 

p 

3.41 

K 

1.53  xlO  3 

q  [m2/s3]a 

lxlO4 

VR  [m2] 

6.7 

oc,  (3 ,  and  K  have  no  dimensions; 


a  Units  according  to  [38], 

For  this  mission,  tracking  by  dead  reckoning  (UUV  navigation  system)  yielded  an 
error  of  476  meters  (see  Figure  70).  The  tracking  algorithm  was  tested  and  the  results 
show  a  more  reliable  predicted  position  for  the  entire  UUV  trajectory  (see  Figure  72). 
When  a  rough  estimate  for  the  distances  was  used,  an  error  of  10.3  meters  was  observed 
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at  the  end  of  the  mission.  While  this  is  a  significant  improvement  over  the  UUV’s 
estimated  position,  it  is  based  on  simple  straight  line  propagation  that  does  not  account 
for  the  refractive  and  multipath  effects  of  the  environment  (Figure  74). 

When  the  more  complete  algorithm  described  in  Chapter  IV  is  used,  the  tracking 
algorithm  produces  an  error  of  8.6  meters  at  the  end  of  the  mission  (as  represented  in 
Figure  73,  there  are  uncertainties  associated  with  the  predictions).  It’s  a  small 
improvement  in  this  case  but  still  indicates  that  correlating  the  measured  impulse 
response  of  the  channel  with  a  ray  trace  prediction  may  provide  an  improvement  in  the 
tracking  accuracy. 


Table  6.  Mission-3  Tracking  Comparison 


Mission-3 

Distances  according  to 
Chapter  IV 

Distances 
roughly  estimated 

Objective  Function  [m4] 

1145 

1742 

Error  at  beginning 
of  the  Mission  [m] 

2.45 

3.10 

Error  at  end 
of  the  Mission  [m] 

8.60 

10.3 

Total  Error  [m] 

11.1 

13.4 
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Figure  74.  Mission-3,  Tracking  Results  when  Processing  the  Distances  Estimated 

According  to  Chapter  IV  and  when  Processing  the  Distances  Roughly  Estimated 
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The  sea  current  predicted  by  the  tracking  algorithm  for  this  mission  is  0.34  knots 
(i.e.,  the  average  of  the  smoothed  predictions)  and  has  a  general  direction  as  represented 
in  Figure  75a.  It  is  worth  noting  the  evolution  of  the  tide  during  the  mission  period  (see 
Figure  75b). 


1400  1600  1800  2000  2200  2400  2600 

x  [m] 


■  M3  12:05  to  13:07  hs 


Figure  75.  Mission-3,  Tracking,  Sea  Current,  and  Tide 


D.  GENERAL  REMARKS  ON  TRACKING  RESULTS 

In  this  section,  a  basic  check  of  the  corrections  provided  by  the  matched  peak  IR 
algorithm  is  made.  Additionally,  a  discussion  about  the  optimum  tuning  parameters  and 
about  the  estimated  sea  currents  is  provided. 

1.  Distance  Corrections 

In  the  beginning  of  this  chapter,  the  term  distance  roughly  estimated  was  defined 
as  being  the  result  of  the  multiplication  of  the  one-way  travel  time  by  the  characteristic 
medium  sound  speed,  corrected  by  the  assets’  depth  (horizontal  distance).  Let  us  define 
the  term  correction  as  the  difference  between  the  distance  roughly  estimated  and  the 
distance  estimated  using  the  approach  defined  in  Chapter  IV.  Table  7  summarizes  the 
amount  of  correction  in  each  of  the  UUV’s  missions. 
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Table  7.  Summary  of  the  Corrections  on  the  UUV’s  Mission 


Number  of 
Measurements  (NM) 

Number  of  Measurements 
with  corrections  >  5  m  (NM5) 

Ratio 

NM5/NM 

Average 
Correction  [m] 

Mission-1 

40 

30 

75.0% 

9.4 

Mission-2 

11 

6 

54.5% 

7.9 

Mission-3 

18 

7 

38.9% 

4.0 

From  Table  7,  it  can  be  seen  in  mission-1  that  75%  of  the  measurements  had  a 
correction  greater  than  5  meters  and  the  average  correction  was  9.4  meters.  Other 
missions  show  smaller  numbers.  The  reason  for  the  difference  between  missions  may  be 
explained  in  the  following  manner. 

When  the  UUV  is  close  to  a  reference  point  and  operating  at  greater  depths,  the 
direct  path  generally  produces  a  strong  first  arrival.  Therefore,  the  algorithm  may  produce 
a  small  correction.  On  the  other  hand,  when  the  UUV  is  far  from  the  reference  point  and 
operating  at  shallow  depths,  the  bottom  bounce  path  dominates  and  a  weak  first  arrival  is 
expected.  Therefore,  the  algorithm  may  produce  a  larger  correction. 

The  ratio  distance  from  the  UUV  to  a  reference  point  ( distance )  over  UUV  depth 
(depth)  can  be  used  as  an  indicator  for  those  situations.  A  smaller  ratio  represents  a 
situation  where  the  UUV  is  closer  to  a  reference  point  and  deep;  therefore,  small  and  less 
frequent  corrections  are  expected.  A  larger  ratio  represents  a  situation  where  the  UUV  is 
farther  from  a  reference  point  and  shallow;  therefore,  larger  and  more  frequent 
corrections  are  expected. 

Table  8  represents  the  average  horizontal  distance  between  the  UUV  and  the 
reference  points  when  a  travel  time  measurement  is  taken,  and  Table  9  represents  the 
average  UUV  depth.  Based  on  these  tables,  the  average  ratio  distance  over  depth  can  be 
calculated. 

The  last  column  of  Table  10  (supported  by  Tables  8  and  9)  shows  the  ratio  of  the 
average  distance  over  depth,  when  the  travel  time  measurements  were  taken  for  each 
UUV  mission.  Mission- 1  has  the  larger  ratio,  followed  by  mission-2  and  mission-3, 
which  means  that  larger  and  more  frequent  corrections  are  expected  in  mission- 1, 
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followed  by  mission-2  and  mission-3.  This  conclusion  is  supported  by  the  data  shown  in 
Table  7. 


Table  8.  Average  Horizontal  Distance  from  the  UUV  to  the  Reference  Points  when 

a  Travel  Time  Measurement  is  Taken 


Avg.  Horizontal 
Distance  from  UUV  to 
R/V  Fulmar  [m]  (DIS) 

Avg.  Horizontal 
Distance  from  UUV  to 
WG  Mako  [m]  (DIS) 

Avg.  Horizontal 
Distance  from  UUV  to 
WG  Tiburon  [m]  (DIS) 

Avg.  Hor. 
Dist.  [m] 
(DIS_A) 

Mission-1 

826 

594 

— 

797 

Mission-2 

610 

840 

252 

497 

Mission-3 

432 

628 

441 

504 

Table  9.  Average  UUV  Depth  when  a  Travel  Time  Measurement  is  Taken 


Avg.  UUV  depth 
when  a 

measurement  is 
taken  to  R/V  Fulmar 
[m]  (DEPTH) 

Avg.  UUV  depth 
when  a 

measurement  is 

taken  to  WG  Mako 
[m]  (DEPTH) 

Avg.  UUV  depth  when 
a  measurement  is 

taken  to 

WG  Tiburon  [m] 
(DEPTH) 

Avg.  UUV  depth 
when  a 

measurement  is 
taken  [m] 
(DEPTH_A) 

Mission- 

1 

39.8 

34.3 

_ 

39 

Mission- 

2 

16.3 

39.9 

33.1 

28.2 

Mission- 

3 

33.1 

49.2 

29.7 

34.6 

Table  10.  Ratio  Average  Horizontal  Distance  to  a  Reference  Point  over  Average 
UUV  Depth  when  a  Travel  Time  Measurement  is  Taken 


Ratio  DIS /DEPTH 
for  R/V  Fulmar 

Ratio  DIS /DEPTH  for 
WG  Mako 

Ratio  DIS  /  DEPTH  for 
WG  Tiburon 

Ratio  DIS_A  / 
DEPTH_A 

Mission- 

1 

20.8 

17.3 

_ 

20.4 

Mission- 

2 

37.4 

21.1 

7.6 

17.6 

Mission- 

3 

13.1 

12.8 

14.8 

14.6 

In  Table  1 1,  the  time  between  the  last  travel  time  measurement  and  the  UUV  GPS 


fix  at  the  end  of  the  missions  is  shown.  As  can  be  seen,  mission-3  has  the  larger  number. 

This  larger  time  difference,  together  with  the  reduced  number  of  measurements,  account 
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for  the  larger  confidence  ellipse  associated  with  the  last  position  prediction,  in 
comparison  with  the  previous  missions. 


Table  11.  Time  between  the  Last  Measurement  and  the  UUV  GPS  Fix 

at  the  End  of  the  Missions 


Time  between  last 
measurement  and 

GPS  fix  [mm:ss] 

Mission-1 

5:43 

Mission-2 

2:21 

Mission-3 

6:37 

2.  The  Tuning  Parameters 

The  tuning  parameters  for  the  UKF-based  tracking  algorithm  can  be  divided  into 
two  groups.  The  first  group  is  composed  of  the  hyper-parameters  a ,  f3 ,  and  k  , 
appearing  exclusively  in  UKF  applications.  The  second  group,  common  to  all  KF 
applications,  is  composed  of  the  matrices  R  and  Q,  representing  the  covariance  of  the 
measurement  and  covariance  of  the  plant  noise  (also  known  as  process  noise), 
respectively. 

As  can  be  seen  in  Chapter  II,  Equation  (2.26),  the  measurements  are  scalar 
numbers,  which  makes  R  a  scalar,  too.  The  square  root  of  the  scalar  R  represents  the 
uncertainty  associated  with  the  measurements  and,  in  our  application,  it  combines  the 
uncertainties  in  the  estimated  horizontal  distances  and  in  the  position  of  the  reference 
points  (in  x  and  y).  The  structure  of  the  covariance  matrix  Q  is  shown  in  Chapter  V, 
Equation  (5.26),  where  q  is  used  for  tuning. 

The  hyper-parameters  a,  (3 ,  and  k  are  scalars  numbers  where  a  and  k 
determine  the  spreading  of  the  sigma  points  around  the  predictions  (mean)  and  f3  is  used 
to  incorporate  prior  information  about  the  non-Gaussian  distribution  of  the  predicted 
states  ( (3  is  set  to  2  to  model  Gaussian  distributions). 
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Table  12  shows  the  optimum  tuning  parameters  for  all  three  UUV  missions.  As 
can  be  seen,  the  parameters  vary  from  mission  to  mission  and,  as  expected,  there  is  no 
one  set  of  tuning  parameters  that  works  universally.  What  draws  our  attention  are  the 
large  variations  between  missions.  Those  variations  cannot  be  attributed  only  to  changes, 
such  as  noise  and  sea  state,  in  the  environment. 

The  problem  solved  by  the  tracking  algorithm  is,  essentially,  a  triangulation 
problem  that  consists  of  range-only  measurements  not  taken  simultaneously.  As  is 
common  in  triangulation  problems,  the  number  of  reference  points,  their  relative 
positions  with  respect  to  the  UUV,  and  the  number  of  measurements  all  influence  the 
tracking  dynamics  and  performance. 

As  seen  in  Figures  60a,  66a,  and  72a,  the  relative  position  between  assets  and  the 
measurement  dynamics  are  very  different  when  the  missions  are  compared.  Those 
differences  may  be  considered  the  major  factor  for  large  variations  in  the  tuning 
parameters,  as  shown  in  Table  12. 

Simulation  shows  that  when  the  measurements  are  consistent  and  in  large 
numbers,  the  influence  of  the  asset’s  relative  position  on  the  tuning  parameters  is 
reduced.  In  this  situation  the  tuning  parameters  for  different  geometries  tends  to  agree, 
which  could  indicate  that  real-time  tracking  is  still  a  possibility.  Although,  as  shown  in 
this  chapter,  with  a  reduced  number  of  measurements,  post  processing  represents  the  best 
choice  to  reconstruct  the  UUV  trajectory. 


Table  12.  Summary  of  the  Tuning  Parameters 


Optimum  Tuning  Parameters 

Mission-1 

Mission-2 

Mission-3 

a 

1.08x10'' 

3.05x10' 

1.42x1  O'2 

fi 

3.82 

8.98xl0'2 

3.41 

K 

2.03  xlO'9 

9.4x1 0"3 

1.53xl0'3 

q  [m2/s3] a 

2.09  xlO'3 

l.llxlO'3 

lxlO'4 

VR  [m2] 

4.65x10s 

10 

6.7 

a,  /?  ,  and  K  have  no  dimensions; a  Units  according  to  [38] 
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3. 


Sea  Current 


The  tracking  errors  at  the  end  of  the  missions,  as  well  as  the  average  predicted  sea 
current,  are  presented  in  Table  13.  As  can  be  seen,  there  is  a  correlation  between  current 
and  the  errors,  where  higher  currents  are  associated  with  higher  errors.  This  problem  is 
common  to  all  underwater  navigation  systems  that  are  not  able  to  estimate  or  measure  the 
sea  current. 

Some  UUVs,  equipped  with  Doppler  velocity  logs,  have  very  successful  results 
measuring  the  sea  current  by  tracking  its  motion  with  respect  to  the  bottom  [53]. 
Unfortunately  this  method  does  not  work  for  deeper  waters  (300  m  and  beyond). 
Research  in  the  area  of  sea  current  estimation  and  measurement  is  growing  very  fast  and 
some  are  leaning  towards  KF  techniques,  as  used  in  this  work  [54],  [55], 


Table  13.  Summary  of  Mission  Errors  and  Average  Sea  Current 


Approximated  Distance 
Traveled  (D)  [m] 

Final  Trajectory 
Error  (E)  [m] 

E/D 

Predicted  Average 
Current  [knots] 

Mission-1 

1400 

1095 

0.78 

0.60 

Mission-2 

1000 

747 

0.75 

0.45 

Mission-3 

770 

476 

0.62 

0.34 

a.  General  Direction  of  the  Sea  Current 

In  an  open  area,  where  the  direction  of  the  water  flow  is  not  restricted  by  any 
barriers,  the  tidal  current  is  rotary  (i.e.,  clockwise  in  the  northern  hemisphere),  flowing 
continuously  throughout  all  points  of  the  compass,  with  speeds  usually  varying  during  the 
tidal  cycle  [56].  The  changes  in  direction  and  speed  may  be  represented  by  arrows  and 
generally  have  an  elliptical  pattern.  On  the  other  hand,  in  rivers  or  straights  the  current 
flows  in  approximately  opposite  directions  with  speed  changing  throughout  the  tidal 
cycle  from  zero  to  a  maximum. 

As  can  be  seen  in  Figure  19,  the  area  of  operation  is  not  exactly  open  and  at  the 
same  time  is  not  extremely  restricted  as  in  a  straight  or  river.  Therefore,  some  influence 
of  the  canyon  (due  its  proximity)  and  of  the  geographical  shape  of  the  bay  may  be 
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expected.  The  tracking  algorithm  predictions  for  sea  current  and  the  tide  information,  for 
all  three  missions,  are  summarized  in  Figures  76  and  77. 

In  Figure  77  it  can  be  seen  that  mission- 1  took  place  on  one  side  of  the  tide  and 
missions  2  and  3  on  the  other  side.  The  change  in  direction  is  clockwise  and  smooth  from 
mission-2  to  mission-3  (as  expected  for  an  open  area),  but  drastic  from  mission-1  to 
mission-2.  This  drastic  change  in  direction  may  be  caused  by  the  change  in  tide  from 
flood  to  ebb  associated  with  some  influence  of  the  canyon  and  the  bay’s  geographical 
shape. 


General  Direction  of  the  Current  -  Origin  Lat:36.79  Lon:-121.89 
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Average  Sea  Current 

Ml :  0.60  knots 
-  M2:  0.45  knots 
M3:  0.34  knots 


_ 

Mission-1 

Mission-2 

Mission-3 

Figure  76.  General  Direction  of  the  Sea  Current  for  the  Three  UUV  Missions 


Date/Tine  (LST/LDT) 


Ml  09:34  to  10:36  hs 
M2  10:49  to  11:52  hs 
M3  12:05  to  13:07  hs 


Figure  77.  Tide  Evolution  during  the  UUV  Mission 
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b.  Initialization  of  the  Tracking  Algorithm  with  Previous  Knowledge  of  the 

Current 

Let  us  consider  a  situation  in  which  a  UUV  has  previous  knowledge  of  the  sea 
current.  We  could  imagine  the  scenario  described  in  Chapter  V,  Section  G,  where  a  fleet 
of  UUVs  is  able  to  transmit  to  the  Command  Center  (via  surface  assets)  their  predictions 
for  the  sea  current,  and  the  Command  Center  (running  the  consensus  current  algorithm) 
broadcasts  back  to  the  UUVs  the  CSC.  As  an  example,  let  us  consider  the  data  from 
mission-2. 

When  mission-2  is  started,  the  UUV  knowing  the  CSC,  initializes  its  tracking 
algorithm  with  this  information.  As  the  first  measurement  did  not  happen  until  the  middle 
of  the  mission,  the  initial  predictions  considering  previous  knowledge  of  the  current  may 
be  improved.  As  can  be  seen  in  Figures  78  and  79,  comparing  the  black  and  the  green 
curves  before  the  first  measurement  (first  “jump”),  the  non-smoothed  tracking  in  green  is 
more  consistent  with  the  final  result  (smoothed  trajectory  in  red)  than  the  black  curve 
where  no  previous  knowledge  of  the  current  was  assumed. 

For  this  particular  case,  due  to  the  geometry  involved  (relative  position  of  the 
UUV  and  the  surface  assets)  and  due  to  the  limited  number  of  measurements,  the  final 
result  was  not  affected  by  the  initial  assumption  for  the  current.  In  a  different  scenario, 
however,  the  initial  knowledge  of  the  current  may  improve  the  overall  performance  of  the 
tracking  algorithm  or,  at  least,  provide  a  faster  convergence  to  the  true  trajectory. 
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08/12/2015  Tests,  Mission-2  -  Origin  Lat:36.79  Lon:-121.89 


Figure  78.  Tracking  Considering  a  Previous  Knowledge  of  the  Sea  Current 
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Figure  79.  Tracking  Considering  a  Previous  Knowledge  of  the  Sea  Current, 

Alternative  View 
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VIII.  CONCLUSION 


Parts  of  this  chapter  were  previously  published  by  EEEE  [1].* 

How  can  a  UUV  be  kept  underwater  for  longer  periods  while  maintaining 
accurate  positioning?  This  is  a  question  that  this  dissertation  intended  to  address.  This 
problem  arises  due  to  the  nature  of  underwater  navigation  where  GPS  signals  (or  similar) 
and  other  navigation  aids  are  not  available. 

Underwater  navigation  typically  relies  on  automatic  dead  reckoning,  gathering 
information  from  several  sensors.  Some  systems  count  on  lower  error,  expensive  inertial 
sensors  and  others  rely  on  simple  magnetic  compasses  and  speed  estimates.  Independent 
of  how  complex  the  system  is,  errors  accumulate  over  time  and  may  lead  to  unacceptable 
position  uncertainties.  Therefore,  a  position  update  from  an  external  reliable  source  is  a 
real  necessity. 

As  pointed  out  in  the  Chapter  I,  by  surfacing,  the  UUV  can  get  a  position  update 
using  its  GPS.  But  this  can  introduce  delays  in  the  UUV’s  mission.  Updating  the  UUV’s 
position  while  it  is  underwater  during  the  course  of  its  mission  may  reduce  the  delays, 
permitting  longer  missions  while  maintaining  position  accuracy. 

State-of-the-art  techniques  estimate  a  UUV’s  position  by  measuring  its  distance, 
and  in  some  systems,  its  bearing  to  a  reference  point  (or  more  than  one)  located  at  the 
surface  or  sea  floor  where,  generally,  the  reference  points — as  well  as  the  UUV — are 
slow  moving  (or  static).  By  measuring  the  travel  time  of  an  acoustic  wave,  a  system  can 
estimate  the  distance  (and  in  more  advanced  systems  the  bearing,  too)  and  the  UUV 
position  can  be  determined  and  updated. 

This  approach  is  used  in  systems  such  as  SBL,  LBL,  USBL,  GIB,  and  some 
hybrid  systems  based  on  the  previous  ideas.  Those  systems  are  reliable  when  the  depth 

*  ©  2016  IEEE.  Personal  use  of  this  material  is  permitted.  Permission  from  IEEE  must  be  obtained  for 
all  other  uses,  in  any  current  or  future  media,  including  reprinting/republishing  this  material  for  advertising 
or  promotional  purposes,  creating  new  collective  works,  for  resale  or  redistribution  to  servers  or  lists,  or 
reuse  of  any  copyrighted  component  of  this  work  in  other  works. 
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being  explored  is  on  the  order  of  or  larger  than  the  horizontal  distances  between  the 
assets. 

As  shown  in  Chapter  IV,  estimation  of  horizontal  distances  in  relatively  shallow 
water  is  a  difficult  task  mainly  due  to  the  multipath  nature  of  the  propagation.  Multiple 
arrivals  reach  the  receiver  at  different  times  by  different  propagation  paths,  making  it 
difficult  to  estimate  accurately  the  travel  time  of  the  acoustic  signals  and,  consequently, 
the  distance.  Some  of  the  systems  described  earlier  are  able  to  account  for  the  refraction 
of  the  sound  waves  when  the  sound  speed  profile  is  available,  but  none  of  them  provides 
a  solution  to  treat  the  errors  caused  by  the  multipath  in  their  distance  estimations. 

Due  to  their  relative  low  cost,  compact  design,  and  reliability,  we  are  relying  on 
battery  powered  DSP-based  acoustic  modems  that  make  use  of  acoustic  communication 
protocols  to  measure  the  acoustic  wave  travel  time  from  the  UUV  to  the  reference  points 
(Chapter  IV,  Section  A).  Chapter  III  presents  the  theoretical  bases  for  wave  travel  time 
measurements  and  the  main  characteristics  of  the  waveforms  used  for  this  task  were 
highlighted. 

Based  on  the  travel  times  measured  by  the  acoustic  modems,  we  developed  an 
algorithm  to  take  into  account  the  multipath  and  the  refraction  of  the  sound  waves  when 
estimating  distances  underwater  (Chapter  IV,  Section  B).  This  algorithm  makes  use  of  a 
ray  tracing  code  to  model  the  environment  and  an  iterative  routine  to  match 
measurements  with  synthetic  predictions.  The  whole  process  is  designed  to  take  just  a 
few  seconds  to  converge  to  a  solution,  making  it  appropriate  for  real-time  applications. 
This  approach  is  a  unique  contribution  of  this  dissertation. 

At  this  point,  using  the  estimated  distances  and  the  coordinates  of  the  reference 
points,  we  may  establish  the  UUV  position.  To  accomplish  this,  it  was  necessary  to 
develop  a  tracking  algorithm  to  fuse  all  available  data  (i.e.,  UUV  navigation  data, 
reference  points’  coordinates,  and  distances  from  the  UUV  to  the  reference  points).  This 
was  done  using  KF  techniques. 

A  model  for  the  UUV  based  on  the  state  space  representation  was  developed  in 
Chapter  II,  taking  into  account  the  two  main  characteristics  of  our  acoustic  modems: 
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inability  to  take  simultaneous  measurements  to  different  reference  points  and  lack  of 
bearing  (direction)  measurement  associated  with  the  travel  time. 

The  aforementioned  modem’s  characteristics  produce  a  non-linear  measurement 
equation  that  was  treated  using  two  different  approaches:  a  linearization  by  Taylor  series 
expansion  as  in  the  EKF,  and  a  statistical  linearization  based  on  what  are  called  “sigma 
points,”  as  in  the  UKF.  To  produce  realistic  trajectories,  mainly  when  the  number  of 
measurements  is  limited  and  during  the  transient  portion  of  the  tracking,  a  smoothing 
algorithm  had  to  be  used,  and  the  equations  are  presented  in  Chapter  V,  Section  F. 

Another  important  part  of  the  tracking  algorithm  is  the  choice  for  the  tuning 
parameters.  In  Chapter  V,  Section  E,  an  algorithm  was  presented  to  select  the  optimum 
tuning  parameters  automatically  for  a  given  set  of  data  (training  data).  This  algorithm 
used  the  UKF  (or  EKF)  to  construct  a  cost  function  and  optimization  techniques  to  find 
the  best  (or  optimum)  set  of  tuning  parameters  that  minimized  (or  maximized)  the  cost 
function. 

To  account  for  the  drift  in  the  UUV  trajectory  caused  by  the  sea  current,  it  was 
modeled  as  a  random  walk  and  was  part  of  the  state  of  the  system.  In  a  scenario  described 
as  a  multi-UUV  centralized  network  (Chapter  V,  Section  G),  the  predictions  for  the  sea 
current  coming  from  the  UUVs  are  concentrated  in  the  command  center,  where  an 
algorithm  processing  this  data  tries  to  establish  a  consensus  current.  When  the  results  of 
the  algorithm  converge,  the  CSC  is  achieved  and  may  be  used  in  several  different  ways  as 
described  in  Chapter  V,  Section  G. 

At  this  point,  one  can  see  that  it  was  necessary  to  put  together  tools  already 
available  as  well  as  to  make  new  developments  to  respond  to  the  short,  but  not  easy, 
question  posed  in  the  beginning  of  this  chapter.  The  first  test  for  this  whole  package  was 
the  processing  of  synthetic  data. 

In  Chapter  VI  a  scenario  in  which  two  UUVs  navigating  in  an  area  where  three 
surface  assets  are  present  was  emulated.  The  surface  assets  were  assumed  to  have  access 
to  GPS  and  were  able  to  share  information  with  the  shore-based  command  center  via 
satellite,  permitting  evaluation  of  the  consensus  current  algorithm. 
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Three  distinct  cases  were  simulated.  The  first  was  the  ideal  case  in  which  all  the 
measurements  were  successful,  and  no  noise  was  present.  The  second  case  considered 
failure  in  the  measurements  during  a  certain  portion  of  the  UUV’s  mission.  In  the  third, 
noise  was  added  to  the  previous  case.  Results  showed  the  importance  of  the  smoothing 
algorithm  to  produce  realistic  trajectories  and  highlighted  that  the  UKF  produced  faster 
convergence  and  better  behavior  in  the  presence  of  noise.  The  consensus  current 
algorithm  proved  converging  in  all  three  cases  with  faster  convergence  and  smoother 
profiles  when  UKF  was  used. 

To  validate  the  developed  algorithms,  real  data  collected  during  the  sea  tests  that 
took  place  in  Monterey  Bay  on  August  2015  was  processed  (Chapter  VII).  In  a  one-week 
sea  test,  two  Liquid  Robotics  Wave  Gliders  (USVs)  and  a  command  ship,  all  equipped 
with  Teledyne-Benthos  acoustic  modems  (ATM-900  series)  and  GPS  connectivity, 
formed  a  network  of  reference  points  for  an  Exocetus  Coastal  Glider  (UUV),  also 
equipped  with  the  same  type  of  acoustic  modem.  On  three  (one-hour  long)  UUV 
missions,  successful  data  was  recorded. 

During  the  course  of  mission- 1,  40  measurements  from  the  UUV  to  the  reference 
points  were  recorded.  The  measurements  were  taken  in  an  inconsistent  manner  and  to 
only  two  of  the  three  reference  points.  After  running  the  automatic  tuning  algorithm,  the 

optimum  value  for  the  parameter  \[r  was  very  high  (see  Table  1).  This  parameter 
combines  the  uncertainties  associated  with  the  distances  between  UUV  and  the  reference 
points,  and  the  uncertainties  in  the  position  of  the  reference  points  (in  x  and  y). 

A  large  ^//?  indicates  a  non-reliable  or,  in  this  case,  inconsistent  measurement. 
The  geometry  (i.e.,  the  relative  position  between  assets),  the  large  period  of  time  with 
measurements  to  only  one  reference  point,  and  measurements  to  only  two  reference 

points  may  be  causing  a  high  ^/R  (see  Figure  59). 

In  mission-2  and  mission-3  fewer  measurements  were  recorded,  but  they  were 
recorded  in  a  more  consistent  way  and  for  all  the  three  surface  assets  (see  Figures  65  and 

71).  That  may  be  the  cause  of  the  reduction  in  the  tuning  parameter  s/~R  (chosen  by  the 

automatic  tuning  algorithm)  to  reasonable  levels  (see  Table  12). 
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As  seen  in  mission- 1,  the  GPS  fixes  for  mission-2  and  mission-3,  that  were 
recorded  as  soon  as  the  UUV  surfaced,  are  located  inside  the  uncertainty  ellipses 
associated  with  the  last  tracking  algorithm  prediction.  It  is  also  worth  noting  that  those 
ellipses  have  realistic  dimensions  (see  Figures  61,  67,  and  73). 

Table  12  showed  the  optimum  tuning  parameters  for  all  three  UUV  missions.  In 
this  table,  unexpected  large  differences  between  missions  can  be  seen.  Those  differences 
may  be  primarily  attributed  to  different  geometries  (i.e.,  the  relative  positioning  between 
assets)  and  a  different  dynamic  in  the  measurements  (i.e.,  the  time  between 
measurements,  total  number  of  measurements,  and  number  of  reference  points  available). 

Simulation  showed  that  when  the  measurements  are  consistent  and  in  large 
numbers,  the  influence  of  the  asset’s  relative  position  on  the  tuning  parameters  is 
reduced.  In  this  situation  the  tuning  parameters  for  different  geometries  tends  to  agree, 
which  could  indicate  that  real-time  tracking  is  still  a  possibility.  Although  with  a  reduced 
number  of  measurements,  post  processing  represents  the  best  choice  to  reconstruct  the 
UUV  trajectory. 

Another  interesting  point  emerged  from  the  predictions  for  the  sea  current.  Those 
predictions  for  all  three  missions  are  depicted  in  Figure  76  (as  the  average  of  the 
smoothed  predictions)  as  well  as  the  evolution  of  the  tide  in  Figure  77.  Note  that  mission- 
1  took  place  on  one  side  of  the  tide  and  missions  2  and  3  on  the  other  side.  The  drastic 
change  in  the  direction  of  the  current  occurred  from  mission- 1  to  mission-2  (and  mission- 
3)  may  be  explained  by  the  change  in  tide  from  flood  to  ebb  associated  with  some 
influence  of  the  canyon  and  the  bay’s  geographical  shape. 

Finally,  to  point  out  the  possible  benefits  of  the  CSC,  data  from  mission-2  was 
used  in  a  scenario  that  simulated  a  broadcast  of  the  CSC  to  all  UUVs  in  a  certain  area. 
Results  showed  that  previous  knowledge  of  the  current  may  improve  the  tracking,  mainly 
in  the  portion  where  few  measurements  are  available. 

Despite  these  promising  results,  improvements  can  be  made  in  several  parts  of  the 
system,  including: 
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-  Measuring  the  beam  patterns  of  the  acoustic  modem’s  transducers  when 
they  are  mounted  in  their  respective  structures  and  including  this 
measurement  in  the  modeling  of  the  system; 

-  Improving  the  corrections  for  the  displacement  between  the  GPS  antenna 
and  the  acoustic  modem  transducer  on  the  Wave  Gliders; 

-  Reducing,  on  the  Wave  Gliders,  the  time  between  GPS  fixes  to  avoid  the 
use  of  interpolation  to  establish  their  position  during  the  measurements; 

-  Accounting  for  the  variability  in  the  sound  speed  profile  measured  by  the 

uuv. 

Continuing  with  a  measurement  system  able  to  provide  range  only,  a  study  of  the 
optimum  tuning  parameters  for  different  geometries,  measurement  dynamics,  ambient 
noise,  and  sea  states  is  an  interesting  topic  for  future  work.  Along  this  line  of  thought, 
different  trajectories  taken  by  a  mobile  reference  point  could  be  explored. 

Although  representing  a  very  interesting  topic,  the  natural  evolution  of  this 
research  points  to  the  use  of  directional  acoustic  modems.  This  equipment,  available  off 
the  shelf,  is  able  to  provide  a  bearing  associated  with  the  travel  time  and  basically  has  the 
same  dimensions  and  weight  as  the  non-directional  version. 

The  use  of  directional  modems  may  reduce  the  number  of  reference  points  to,  in 
theory,  only  one.  It  will  eliminate  the  influence  of  different  geometries  and  measurement 
dynamics  in  the  optimum  tuning  parameters,  although  ambient  noise  and  sea  state  may 
still  affect  them. 
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APPENDIX  A.  MATCHED  FILTER  OUTPUT  FOR  LFM  AND  HFM 
PULSES  IN  THE  PRESENCE  OF  DOPPLER  EFFECT 


The  purpose  of  this  appendix  is  to  demonstrate  how  the  relative  motion  between 
source  and  target  can  affect  the  matched  filter  output  for  LFM  and  HFM  pulses. 

A.  LFM  PULSES 

Let  us  start  with  the  evaluation  of  the  Fourier  transform  of  the  rectangular- 
envelope  LFM  pulse. 


Fourier  Transform 


Let  us  start  with  the  rectangular-envelope  LFM  pulse  in  complex  notation 


x(t)  =  A  rect 


f  *  \ 


J 


i{lnfct+Dpt 2) 


(A.l) 


TC 

where  D  =  ±  —  B ,  and  rect 


't' 


V  j 


=  1  for  t  <  T/2 .  The  Fourier  transform  of  x(t)  is 


X  (/)  =  J  A  rect 


J 


ej(2,fct+D/)  e_J2xfidt  ' 


(A.2) 


Rearranging  the  terms  in  Equation  (A.2),  we  obtain 


» *o 


X(f)=  j  A  rect 


f  *  \ 
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j(2x(fc-f)t+Dpt2 


dt . 


(A.3) 


A(t) 


To  evaluate  Equation  (A.3)  the  method  of  stationary  phase  [14],  [57],  and  [58]  is 
used.  The  method  of  stationary  phase  is  useful  to  make  approximations  for  integrals  with 
highly  oscillatory  integrands.  The  stationary  phase  point  is  defined  such  that  the  first  time 
derivative  of  the  phase  in  Equation  (A.3)  is  set  to  zero,  ^'(t0)  =  0.  Then  the 
approximation  of  the  integral  of  Equation  (A.3)  is  [14] 


x(f)' 


-2  n 


The  integrand  phase  and  its  derivatives  are 


AM 


(A.4) 
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<j>(t)  =  2x(fr-f)t  +  D/ 
f(t)  =  2n(fc-f)  +  2Dpt 

r(t)  =  2D„. 

The  stationary  phase  point  is  defined  by 

2x(fc-f)  +  2Dpt0  =0 

U P 

Using  the  definition  of  D  and  defining  k  =  B/T ,  we  can  write 


to=j(f~fc)- 


(A.5) 


(A.6) 


(A.7) 


Using  the  definitions  of  A(t)  and  (J){t)  from  Equation  (A.3)  and  substituting  Equation 
(A.7)  into  Equation  (A.4),  we  obtain 

X  (/)  *  eJ *  rect^^j  e'Jt(/”/c)  1 .  (A.8) 

Equation  (A.8)  is  the  approximate  Fourier  transform  of  the  rectangular-envelope  LFM 
pulse  by  means  of  the  stationary  phase  method.  This  result  will  be  used  next,  in  the 
evaluation  of  the  matched  filter  output. 


2.  Matched  Filter  Output  in  the  Presence  of  Doppler  Effect 

The  following  development  is  based  on  the  methodology  described  by  Minkoff 
[17].  Re-stating  Chapter  III,  Equation  (3.17),  the  Doppler  transformed  received  signal  is 

x\t)  =  x{at-T0),  (A.9) 

/  \  2v 

where  x[t)  is  the  transmitted  signal,  z0  is  the  two-way  travel  time,  and  ««1 - . 

c 

Worth  noting  that  vr  is  the  relative  radial  speed  between  source  and  target  as  described  in 
Chapter  III,  section  B.l. 

The  output  of  a  filter  matched  to  the  transmitted  signal  x(t)  is 

y(t)  =  x(orr-ro)*x*(-r),  (A.10) 

where,  in  Equation  (3.6),  t0  was  set  to  zero  and  K  to  1. 

Defining 


140 


F(/)  =  F{x(«r-ro)*x*(-r)} 

.  .  1  (  f\  -j2x—r  *.  . 
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(A.  12) 


Using  Equations  (A.8)  and  (A.  12),  the  matched  filter  output  for  the  rectangular- 
envelope  LFM  pulse  is 

.  .  A2  T  r  (  f  la-  f  \  (  f-f') I  -j**L*o  -y{f[(//«-/c)2-(/-/c)2l}  * 

y(t)  = - f  rect  - -  rect  - -  e  “  e  ^  (A.13) 

V  ’  a  Bi  {  B  J  {  B  ) 

where  a  is  always  positive. 

The  exponential  that  contains  the  quadratic  terms  may  be  simplified  as 


-j\j  (//«-/cf-(/-/c)2 
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The  term  that  contains  a 2  in  (A.  14)  may  be  simplified  as 

2  i-(i -Vf  V-4 iyf 
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Considering  that  c»vr  the  quadratic  terms  in  (A.15)  may  be  ignored  and  the 


term  1-  y  « 1 .  Rewriting  (A.15)  yields 


1-a  _  4vr 

Appling  the  same  reasoning  in  - — — ,  we  obtain 

a 


(A.  16) 
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1  -  a  ~  2vr 
a  c 


(A.  17) 


Substituting  (A.  16)  and  (A.  17)  into  (A.  14)  we  obtain 
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Now,  for  further  simplification  in  Equation  (A.25),  let  us  assume  that  fd 

4  f  2 

and  separating  the  constant  termT-Ai^ ,  the  exponential  becomes 
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Substituting  (A.  19)  into  Equation  (A.  13)  yields 
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Rearranging  the  terms  and  rewriting  Equation  (A.20),  we  obtain 


| —  /  Ji\vr  fc: 

y(>)=A2JjA  ‘ 


-rect 


/-/c 

V  ^  y 


71 4v, 


c  k 


1  ' 

—  rect 
a 


f/a-fc)  j2*a\T~T° 


V  5  y 


Equation  (A.21)  may  be  written  as 

y(t)  =  A 


A/) 

(A.21) 
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where  a(t)  =  F  1|A(/)|  and  £>(?)  =  F  1  {#(/)}. 


Referring  to  Equations  (A.l)  and  (A.8) 
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Therefore 
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b(t)  =  |  rect  — ■ 
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Defining  the  following  variables 
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Substituting  (A.26)  into  (A.25)  and  rearranging  the  terms  yields 

Z>(f)  =  e  U  j  j  reef  4-  e  U  j  df'. 


(A.25) 


(A.26) 


(A.27) 


Solving  Equation  (A.27),  using  the  Euler  identity  ei&—e  jd  =  2jsmO ,  and  using 
,  .  sin(^-x) 


dnc(x)  = 
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u 

Substituting  Equations  (A.24)  and  (A.28)  into  (A. 22)  yields 


(A.28) 
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(A.29) 

Taking  the  absolute  value  of  Equation  (A.29),  we  obtain 

\y(t)\  =  A2\fr^B  rect  —  e^2rifcl+Dp '  ^  *  e]2,cfcat  sine  B  —  -  z  +  at 

1  W|  Ir'J  Ik 


(A.30) 


Expanding  Equation  (A.30)  in  the  convolution  integral  format  results  in 
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Adding  and  subtracting  at  in  the  first  exponential  and  rearranging  the  exponential  terms, 
we  obtain 
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(A.32) 

Rewriting  Equation  (A.32)  using  the  convolution  representation,  we  arrive  at 
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Recognizing  that  the  exponential  in  Equation  (A.33)  will  vary  little  over  time  T' 
(t  '  <sc  T,  fd  <sc  fc ,  and  Dp '  <K  Dp  ) ,  we  can  write 
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It  is  commonly  accepted  that  the  width  of  the  convolution  of  two  functions  is 
nominally  the  sum  of  the  widths  of  the  two  functions.  Defining  the  width  of  the  rect 
function  as  AR ,  the  width  of  the  sine  function  as  As ,  and  the  width  of  matched  filter 

output  as  AT  : 

AT~AR+AS.  (A.35) 

Since  the  first  zero  in  the  sine  function  occurs  when  its  argument  goes  to  unity 
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T.  fd 


The  term  —  — represents  a  partial  temporal  offset  of  the  main  lobe  in  the  matched 
a  ka 

filter  output.  Therefore,  the  width  due  to  the  sine  function  is  1/ aB . 


The  width  of  the  rect  function  is  T ' ,  using  Equation  (A.24)  yields 
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Substituting  Equation  (A.37)  and  (A.38)  into  Equation  (A.35),  we  obtain 
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From  Equation  (A.39)  it  can  be  seen  that  in  a  stationary  source  and  target 
scenario,  the  nominal  half-width  of  the  main  lobe  in  the  matched  filter  output  is  1 /B .  Due 
to  the  Doppler  effect,  there  is  a  widening  in  the  main  lobe  that  will  be  accompanied  by  a 
reduction  in  amplitude  (due  to  conservation  of  energy),  as  can  be  seen  in  Chapter  III, 
Figure  9.  Additionally,  there  is  a  temporal  offset  in  the  matched  filter  output’s  peak  that, 
for  some  applications,  must  be  taken  into  account. 


Depending  on  how  severe  the  widening  effect  is,  it  may  be  necessary  to  use  a 
“bank”  of  filters  matched  to  different  target  speeds  to  permit  target  detection.  According 
to  Equation  (A.39)  the  widening  in  the  matched  filter  output  is  given  by  the  factor 
l-i- (4 Vr/C)TB,  and  the  range  resolution  is  degraded  by  the  same  factor.  Thus,  if 

(4 vr/c)ra  <<c  1  the  widening  and,  consequently,  the  amplitude  reduction  in  the  matched 
filter  output  will  not  be  extensive.  Therefore,  it  may  be  written  as 

TB  «  —  .  (A.40) 

4vr 

Equation  (A.40)  states  that  if  the  signal’s  time  bandwidth  product  is  small  in 
comparison  with  the  ratio  of  medium  sound  speed  over  the  maximum  source-target  radial 
speed,  the  widening  effect  in  the  matched  filter  output  can  be  neglected. 

At  this  point  we  may  go  back  and  revisit  Equation  (A.20): 
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Based  on  the  previous  derivation  and  analyzing  the  exponential  term,  we  may 
conclude  that: 


a)  The  first  term  (quadratic),  also  called  the  dispersive  term,  is  responsible 
for  the  widening  in  the  matched  filter  output; 

b)  The  second  term  is  responsible  for  the  temporal  offset  in  the  matched  filter 
output. 

If  we  use  a  waveform  such  that  the  dispersive  term  is  not  present  in  the  matched 
filter  output,  the  widening  effect  as  seen  in  the  LFM  pulse  is  not  expected.  With  this  in 
mind,  next  section  presents  the  analysis  for  the  rectangular-envelope  HFM  pulse. 


B.  HFM  PULSES 

Let  us  start  with  the  evaluation  of  the  Fourier  transform  of  the  rectangular- 
envelope  HFM  pulse. 


1. 


Fourier  Transform 


The  rectangular-envelope  HFM  pulse  may  be  represented  as  [21] 
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where  rect 
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=  1  for  t  <T/ 2;  fa  and  fend  are,  respectively,  the  starting  and  the 
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ending  frequency;  and  T  is  the  duration  of  the  pulse. 
In  complex  notation  Equation  (3.22)  becomes 
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The  Fourier  transform  of  x(t)  is 
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Rearranging  the  terms  in  Equation  (A.43),  we  obtain 
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Applying  the  method  of  stationary  phase  gives  us 
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The  stationary  phase  point  is 
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Using  the  definitions  of  A(?J  and  (J)(t)  from  Equation  (A.44)  and  substituting 


Equation  (A.47)  into  Equation  (A.45),  we  obtain 
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Equation  (A.48)  is  the  approximate  Fourier  transform  of  the  rectangular-envelope  HFM 
pulse  by  means  of  the  stationary  phase  method.  This  result  will  be  used  next,  in  the 
evaluation  of  the  matched  filter  output. 
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2.  Matched  Filter  Output  in  the  Presence  of  Doppler  Effect 

The  following  development  is  based  on  the  methodology  described  by  Minkoff 
[17].  Using  Equations  (A.48)  and  (A.  12),  the  matched  filter  output  for  the  rectangular- 
envelope  HFM  pulse  is 
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Making  the  following  approximations  yields 
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Using  (A.50)  and  rearranging  the  terms,  Equation  (A.49)  becomes 
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Because  the  phase  of  the  exponential  in  (A.51)  is  linear  in  /  ,  the  dispersive  term  is  not 
present.  Therefore,  the  widening  effect  in  the  matched  filter  output  is  not  expected. 


The  function  red 


Vf-Vfo 

.  B/  ■  f of  end 
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~\2 


has  amplitude  one  and  width  defined  as  follows 
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+1<I< 
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1  „  1 


1  B 

-  +  - 


</< 


1 


B 


fo  ^f of  end 

fo 


1  + 


B 


</< 


fo  ^f of  end 

fo 


1- 


B 
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Applying  the  limits  of  the  integral  in  Equation  (A.51)  yields 
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2  l~B/2fend 
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(A.52) 


(A.53) 
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From  Equation  (A.53)  we  want  to  determine: 

a)  The  time  instant  when  the  peak  in  the  matched  filter  output  occurs; 

b)  The  width  of  the  main  lobe  assuming  that  it  has  a  form  of  sine  (x). 

Following  [17],  it  is  not  necessary  to  evaluate  explicitly  Equation  (A.53)  to  answer  (a) 
and  (b). 


For  (a),  considering  that  the  maximum  in  Equation  (A.53)  occurs  when  the 
exponential  goes  to  unity  (which  can  be  thought  of  as  the  stationary  phase  argument),  it 
may  be  written  as 

t  (1  -a 

tP=-~- - 

a  a 

Note  that  (A.54)  is  constant  depending  only  on  the  relative  radial  speed  between  source 
and  target  (factor  a ). 


)  1 


Kfo 


(A.54) 


Assuming  that  l//2  is  slowly  varying,  for  (b)  the  zeros  in  the  sine  (x)  function 
will  be,  in  an  approximation,  considered  occurring  at  an  integral  number  of  cycles, 
leading  to 


where  n  -  ±  1, 2, 3 . . . 


Li 

1  “ 

r  */Jj 

fO(t)  =  n , 


(A.55) 


The  limits  of  the  integral  in  Equation  (A.53)  (to  replace/)  are 

f, _ /,  g(/,//<w) 

l-B/2/w  l  +  S/2/.w  1  ~(B/2f,J' 

Equation  (A.56)  may  be  approximated  as 

B 


Substituting  (A.57)  into  Equation  (A.55)  yields 


B 

l-(2?/2  JZf 


(A.56) 


(A.57) 


(A.58) 
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For  n  =  1  and  rearranging  the  terms  in  Equation  (A.58),  we  can  write 


1 

+  — 

'i-i 

r  B  ' 

2 

L  Kfo  \ 

B 

4 

V  ^ ?nd  ) 

The  first  term  in  Equation  (A.59)  represents  the  center  of  the  main  lobe  in  the 
matched  filter  output  and  the  second  represents  the  half- width.  Since  B/  fend  <  1  always, 

the  width  of  the  main  lobe  in  the  matched  filter  output  is  slightly  different  from  1/ B  and 

independent  of  (4 vr/c)TB.  Thus,  the  widening  and,  consequently,  the  reduction  in 

amplitude  in  the  matched  filter  output  will  not  be  as  severe  as  in  the  LFM  case  (see 
Chapter  III,  Figure  10). 

From  Chapter  III,  Figure  10,  it  can  be  seen  that  with  the  use  of  an  HFM  pulse  the 
filter  is  “always”  matched  to  the  target  speeds,  as  evidenced  by  the  small  decrease  in 
amplitude  as  Doppler  errors  increase,  avoiding  the  use  of  a  “bank”  of  filters.  From 
Chapter  III,  Figure  9,  for  an  LFM  pulse  there  is  a  large  decrease  in  amplitude  as  Doppler 
errors  increase.  In  both  cases  there  is  a  temporal  offset  in  the  matched  filter  output’s  peak 
due  to  the  relative  motion  between  source  and  target. 


In  cases  where  the  relative  radial  speed  (or  Doppler  shift)  is  known,  with  the  use 
of  HFM  pulses,  the  time  of  arrival  may  be  predicted  using  Equation  (A.54)  as  in  the 
following 


(1  -a)  1 

C,*- - Z— r  +  tp ,  (A.60) 

a  KJo 

where  a  « l-2vr/c ,  and  tp  is  the  matched  filter  output  peak  time.  The  temporal  offset  in 
the  matched  filter  output  is 

(l-tf)  1 


t,  =  T  ~t 

d  o  p 


(A.61) 


«  Kf„ 

Therefore,  when  the  Doppler  shift  is  known,  Equation  (A.61)  permits  a  direct 
compensation  for  the  time  offset  in  the  matched  filter.  Unfortunately  there  is  no  closed 
form  expression  for  LFM  pulses.  Continuing  this  discussion,  let  us  analyze  a  simplified 
argument  regarding  the  widening  in  the  matched  filter  output  due  to  Doppler  effect. 
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c. 


DOPPLER-INVARIANT  WAVEFORMS 


Some  authors,  as  in  [21],  [59],  [60],  and  [61],  make  the  argument  (directly  or 
indirectly)  that  if  the  instantaneous  frequency  of  the  transmitted  signal  is  merely  a  time 
delayed  version  of  the  instantaneous  frequency  of  the  Doppler  transformed  received 
signal,  the  widening  and  the  amplitude  reduction  in  the  matched  filter  output  will  not  be 
present,  and  they  label  a  waveform  with  such  characteristics  a  “Doppler-invariant” 
waveform.  Therefore,  for  a  “Doppler-invariant”  waveform  the  following  equation  shall 
be  satisfied 

(A.62) 

where  td  is  a  constant,  f  is  the  instantaneous  frequency  of  the  transmitted  signal,  and 

firec  is  the  instantaneous  frequency  of  the  Doppler  transformed  received  signal.  The  basis 

of  the  preceding  argument  is  that  if  Equation  (A.62)  is  satisfied,  the  quadratic  term  in  the 
matched  filter  output  will  not  be  present.  Thus,  it  does  not  lead  to  the  widening  effects  as 
previously  demonstrated. 


Let  us  start  with  the  LFM  pulse.  From  Equations  (A.  10)  and  (3.19),  the  Doppler 
transformed  received  signal  is 


x'(t)  =  a(at) cos  2 7rfcat+Dp(at)2-T0 


(A.63) 


for  which  the  instantaneous  frequency  is 

firec  (0  =  ^^[27rfcClt  +  DP  ~  T° 


(AM) 


=  fca+- 


cc2D„ 


n 


According  to  Equation  (3.20),  the  instantaneous  frequency  of  the  transmitted  signal  is 


1 


f,(t)  =  f,+-Drt. 


n 


Substituting  Equations  (A.64)  and  (A.65)  into  Equation  (A.62),  we  obtain 


f  ~  2 


*fc 


aAD, 


■(ar-l)  + 


fi-A) 

a  ) 


t . 


(A.65) 


(A.66) 
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Because  td  is  a  function  of  time,  the  LFM  pulse  does  not  satisfy  Equation  (A.62).  Thus, 

the  widening  effect  in  the  matched  filter  output  may  be  expected.  Note  that  there  is  no 
closed  form  expression  to  calculate  the  time  offset  in  the  matched  filter  output’s  peak  for 
an  LFM  signal. 


Moving  to  the  HFM  pulse,  from  Equations  (A.  10)  and  (3.22),  the  Doppler 
transformed  received  signal  is 


2n 

x(r)  =  a  (at)  cos  — ln(l  +  Kfoat)-r0  . 

1C 


(A.67) 


for  which  the  instantaneous  frequency  is 


fo^_ 

-  +  Kf,,at 


(A.68) 


According  to  Equation  (3.23)  the  instantaneous  frequency  of  the  transmitted  signal  is 


(A.69) 


Substituting  Equations  (A.68)  and  (A.69)  into  Equation  (A.62)  yields 

(1  -a)  1 

a  Kf o' 


(A.70) 


In  Equation  (A.70)  td  is  constant  satisfying  the  criteria  for  a  “Doppler-invariant” 

waveform.  Therefore,  the  widening  effect  in  the  matched  filter  output  is  not  expected.  It 
is  worth  noting  that  Equation  (A.70)  agrees  with  Equation  (A.61). 
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APPENDIX  B.  OBSERVABILITY  OF  THE  DYNAMIC  MODEL 


In  this  appendix  we  address  the  issue  of  observability  of  the  overall  system.  In 
particular  we  want  to  make  sure  that  the  dynamic  system  defined  in  Equation  (2.17)  with 
the  observations,  Equation  (2.25),  has  suffient  information  so  that  both  position  and 
current  can  be  estimated. 


Although  the  system  is  clearly  nonlinear,  it  can  be  framed  as  a  linear  time 
invariant  system  with  nonlinear  memoryless  observations.  As  a  consequence  we  can  use 
the  standard  observability  test  using  the  observability  matrix. 

Basic  Fact: 


where  x< 


j  ,Vx  I 


(B.l) 


It  is  well  known  that  given  the  state  space  representation  of  a  dynamic  system 

x  =  Ax  +  Bu 
y  =  Cx 

the  system  is  “observable”,  i.e.,  the  state  x(7)  can  be  estimated  from  u(?)  and  y  (t), 
provided  that  the  following  matrix  (observability  matrix) 


C 

CA 

CA2 


CA 


N-l 


is  full  rank. 


(B.2) 


In  our  case  we  have  the  following: 

-  x  represents  the  position  of  the  UUV,  where  x  e  M2xl ; 

-  c  represents  the  velocity  of  the  sea  current,  where  c  e  M2xl . 

Then  we  can  write 

x 
c 

a)  UUV  and  three  reference  points: 

Suppose  we  have  three  reference  points  PA  ,  PB ,  and  Pc ,  as  depicted  in  Figure  80, 


0 

I 

X 

0 

0 

C 

+  v . 


(B.3) 


we  may  write 
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rA  =|*-PA|. 

rB  =|X-Pb|’ 

r,  =  |x  -  PJ . 


Then  notice  that 


r,2  —  r 2  =|x-PA  T  —  |x  —  P, 


A  B 


B  |  ’ 


as 


a|  —  |b|  =(a+b)  (a-b)  ,  Equation  (B.5)  becomes 

'■i-'i5=(2x-PA-P„)r(PB-PA). 

Similarly 

rl~rc  =  (2x-PA  -Pc)t(Pc-Pa). 
We  can  write  (B.6)  as 


^  ~rl)  =  (Pb  -Pa  f  x-i(PA  +PB)r  (P,  -PA) . 


2  v  ^  d  )  v  **  a  /  2 

Rearranging  Equation  (B.8)  we  obtain 


similarly 


y  ab = \  (’■I  -  ri )  ■ +  i(|p,  r  -  |P*  f  J  =  (P,  -  pA  f x , 
y*c  rl )  +  5(|Pc|2  -|P*  f  f  =  (Pc  -  Pa  f  X  ■ 


2V  L  >  2 
Then  the  overall  state  space  representation  becomes 
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(B.4) 


(B.5) 

(B.6) 

(B.7) 

(B.8) 

(B.9) 

(B.10) 


CB.il> 


To  construct  the  observability  matrix,  the  matrix  A  in  Equation  (B.2)  may  be 
written  as 

“0  I 
0  0 


A  = 


(B.12) 


and  the  matrix  C  as 


C  = 


(P.-Pa)  » 

(Pc-P»)r  # 


(B.13) 
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Using  Equations  (B.2),  (B.12),  and  (B.13)  the  observability  matrix  may  be  written 


as 


C 

CA 


(p,-p*)r  0 

(Pc-P»)r  # 

»  (Pb-P./ 

»  (Pc-P»f 


(B.14) 


As  long  PB-PA  and  PC-PA  are  linearly  independent,  i.e.,  det(S)^0,  the 
system  is  observable. 


X 


Figure  80.  UUV  and  Three  Reference  Points 
b)  UUV  and  two  reference  points: 

As  observability  is  independent  on  the  reference  frame,  let  us  consider  a  scenario 
where  the  origin  of  our  coordinate  system  is  the  reference  point  PA  and  the  line  between 

PA  and  the  reference  point  PB  represents  the  x  axis  as  depicted  in  Figure  81. 
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y  t 

yu- 


PB 

0 - ► 

x 

Figure  81.  UUV  and  Two  Reference  Points 


In  this  case  we  can  write 


and 


p*=[°  of 
p„=[*»  of. 


From  Figure  82  we  can  write 


2  2  2 
rA  =xu  +  yU’ 


i=(xu  ~xh)  +y2u 
=  x2u~^XuXB+x2B  +  yt- 


From  (B.  16)  and  (B.17) 


v  —  v  +  x  —  x 
'  A  'b  ^  ^B  ^~^iiB  • 


From  Figure  8 1  we  can  write 


x„ 


cos  0  =  — , 


substituting  (B.18)  into  (B.19)  we  may  write 

r2_r2+JC2 

cos  6  =  — — - —  =>  #  =  acos 


2rAXB 


2  .  2  ^ 

'5  + 
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lrAXB  J 


(B.15) 

(B.16) 

(B.17) 

(B.18) 

(B.19) 

(B.20) 


From  Figure  8 1  we  may  write 


=tg(6»)x„  ^0  =  tg(6»)x„-ytt.  (B.21) 

From  Equations  (B.18),  (B.20),  and  (B.21)  the  two  observations  may  be  written  as 

$AB 


VA  rB  XB  
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0  =tg(#)xu-y«- 


(B.22) 


156 


Then  the  overall  state  space  representation  becomes 
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X 
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1  0  0 
tg#  -1  0 
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From  (B.2)  and  (B.23)  the  observability  matrix  may  be  written  as 


C 

CA 


10  0  0 

tg  0  -10  0 

0  0  10 

0  0  tg#  -1 


As  det(5)  ^  0 ,  the  system  is  observable. 


(B.23) 


(B.24) 
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